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DEFINITION OF SYMBOLS 


Symbol 

a 

A 

C P 

f 


m test 

p 


p c 

p 


Q 

R 

s 


t 

T 

u 


7 

P 




Definition 

nondimens ional speed of sound a*/a* 

cross-sectional area 

specific heat at constant pressure 

f * L 0 

nondimens ional quantity — — where f ,f is the sum of body 

a * 2 

o 

and dissipative forces per unit mass 
total length of facility (with dimensions) 

Mach number of shock 

Mach number in test section during steady state duration 

nondimens ional static pressure p*/p* 

initial pressure ratio across diaphragm 

right running characteristic variable 

left running characteristic variable 

gas constant 

nondimens ional specific entropy s/Cp( 7 -l) 
nondimens ional time aot' v 7L 0 
nondimens ional temperature T' v /Tq 
nondimens ional flow velocity u'Va^ 
nondimens ional shock velocity w^/a^ 
ratio of specific heats 
nond imens ional dens ity 

L a'' 
o o T 

nond imens ional quantity ; where \|/' v represents the 

7 a p* 

mass flow removed through the walls per unit length 


v 


Subscripts 


DEFINITION OF SYMBOLS (Continued) 
Definition 


o reference conditions which for the case at hand were the 

conditions on the right side of the diaphragm before rupture 

S shock conditions 

C.S. contact surface conditions 


Superscripts 

* 


d imens ional quantities 
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APPLICATION OF THE CHARACTERISTIC METHOD IN CALCULATING THE 
TIME DEPENDENT, ONE -DIMENSIONAL, COMPRESSIBLE FLOW 
IN A TUBE WIND TUNNEL 

SUMMARY 


This report presents in detail a calculation method for one-dimen- 
sional, time-dependent flow through a pressure- tube wind tunnel. The 
computational procedure involves the method of characteristics for the 
one -dimensional unsteady flow of a perfect gas and accounts for area 
changes, shock formations, and intersection of discontinuities in the 
flow field. 

A computer program written in Fortran language was constructed for 
the CDC 3200 digital computer and is presented along with a detailed 
description of the preparation of input for the program. 

Calculated results of the program are presented for a pressure-tube 
wind tunnel now under construction at MSFC. 


I. INTRODUCTION 


Ludwieg [7] first proposed the principle of the pressure tube wind 
tunnel and supervised the construction of the first facility of this type 
at the Aerodynamische Versuchsanstalt in Gottingen. Some experiments 
were reported in reference 8 of the stagnation pressure loss, and some 
schlieren pictures were made of the supersonic jet at the outlet of the 
nozzle. Calculations and measurements were given that show the extent 
of stagnation pressure losses in the test section with increasing bound- 
ary layer thickness. It was concluded that for large tube diameters 
(and therefore large Reynolds numbers) the limit length (stagnation 
pressure loss less than 1 percent due to the boundary layer) of tne tube 
is approximately 100 tube diameters. 

A supersonic pressure tube wind tunnel was constructed at the Royal 
Armament Research and Development Establishment in England during 1957, 
and some measurements from this wind tunnel are reported in reference 4. 
Reference 3 reported further measurements of the pressure at the nozzle 
end of the tube and static and pitot pressures in the working section of 
the nozzle. The diaphragm was located at the nozzle exit. Pictures of 



the flow around a model were made by a high speed camera. An approximate 
calculation procedure for one-dimensional flow reported in reference 3 
was used to obtain analytical results for running times and static pres- 
sures in the test section. Their analytical procedure assumed that the 
expansion fan emanated from the nozzle throat rather than from the 
diaphragm location, and thus neglected the unsteady expansion from the 
diaphragm through the nozzle throat. The effects of cross-sectional area 
changes were accounted for by steady state assumptions. An approximate 
method for calculating the boundary layer growth along the tube is also 
presented. 

Reference 6, which is a good discussion on the solution of hyperbolic 
differential equations that describe one-dimensional, non-steady, com- 
pressive, multi- isentropic flow, presented three methods for solving the 
differential equations, along with the advantages and disadvantages assoc- 
iated with each method. Some hand calculations for an example of isen- 
tropic flow were made, but no extensive application of the procedures was 
performed. 

Bull [1, 2] gives some measurements and calculations of the starting 
processes on an intermittent supersonic wind tunnel that were made at the 
Institute of Aerophysics at the University of Toronto. This tunnel con- 
sisted of a vacuum reservoir at the end of a Laval nozzle with a cello- 
phane diaphragm located either upstream or downstream of the throat. 

Some calculations were made for one case assuming the flow to be time- 
dependent and one-dimensional , and a wave diagram was presented for the 
early phase of the flow from these calculations. 

Rudinger [9] presents an excellent text on calculation procedures 
for solving the partial differential equations of non-steady, one-dimen- 
sional flow of compressible fluids through a duct. Because of the con- 
sistent and straightforward presentation of the solution techniques 
described in this book, many of the calculation procedures described in 
the present report were borrowed from this source. 

Dahm [5] discussed a proposed Ludwieg-tube type of facility at the 
Marshall Space Flight Center for aerodynamic testing, at or near full 
scale Reynolds number, of a Saturn V rocket. This proposal generated 
an interest at the Center in developing an analytical capability for 
calculating start times and other properties of the flow pertinent to 
this facility, thus leading to the material presented here. Some experi- 
mental results of the starting characteristics for a small-scale pilot 
model of a blowdown wind tunnel that was tested at MSFC are presented 
in reference 11. 
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This paper presents a numerical procedure for solving the partial 
differential equations that describe the flow in a pressure- tube wind 
tunnel by a method generally referred to as the f, method of character- 
istics. " The mathematical model of the flow is assumed to be time- 
dependent, one-dimensional, mul ti-isentropic, and compressible. A 
computer program in Fortran IV language was formulated, and wave dia- 
grams and other results of the flow that were calculated from this 
program are given. 

The tunnel considered here (see Figure 1) consists of a long tube 
that serves as the storage reservoir, a convergent-divergent nozzle and 
test section, and an outlet into either the atmosphere or an emptying 
reservoir. The tube, closed at one end, is divided downstream of the 
test section by a diaphragm. When the diaphragm is ruptured, the high 
pressure gas on the left side of the diaphragm expands and compresses 
the gas on the low pressure side, thus creating a shock wave, contact 
surface, and an expansion fan (see Section IV)- The contact surface is 
defined as the gas particles that were initially in contact with the 
diaphragm surface. The expansion fan, bounded on the left by what is 
termed here as the head wave and on the right by the tail wave, is 
mathematically the family of left-running characteristics. The solu- 
tion consists of tracing the left- and right-running characteristic 
curves in the (x,t)-plane after the flow starts at time = 0. 

The starting time is defined as the time when the flow properties 
in the test section become constant. The end test time is defined as 
the time that it takes the headwave to travel along the tube, reflect from 
the closed end of the tube, and reach the test section. Thus, the use- 
ful run time is the difference between the end test time and the start 
time. It is fairly evident that run times can be increased by increas- 
ing the length of the tube. Ludwieg [8] pointed out that to keep 
boundary layer effects negligible the tube length should not exceed 
approximately 100 tube diameters. Boundary layer effects can be 
approximately accounted for in the method presented in this report. 


II. THE DIFFERENTIAL EQUATIONS 

The main purpose of this report is to present the practical numeri- 
* cal procedures for solving the time-dependent equations for one-dimen- 

sional flow through a wind tunnel. Therefore, the fundamental differential 
equations which are derived in numerous references will be listed here 
only in the final form (see reference 9): 


3 



Continuity Equation 


+ SifiVAl + = 0 


Momentum Equation 


Du* 

Dt* 



&* . . jLa£ + 

ox* p* dx* 


f*. 


( 1 ) 


( 2 ) 


Equation of State 


p 

it £_ _ 

=? — *y *-r = 


7 - 7RT . 


Integrated Form of the First Law of Thermodynamics 


3 * - s* = c In — c ^ In ^ . 

o P T * P 7 P* 
o 


( 3 ) 


( 4 ) 


The subscript o indicates some state from which entropy changes are 
measured. The star superscript refers to dimensional quantities. For 
a definition of the symbols, the reader should refer to the Definition 
of Symbols. The underlying assumptions in the derivation of the above 
equations are: 

•j 1 - 

(1) All quantities depend on the time t and a single 
coordinate x*. 

(2) There is only one velocity component u" and that is 
in the x*-direction. 

(3) The gas follows the ideal gas laws, and the values 
of the specific heat are constant. 

(4) All body and dissipative forces are lumped into a 
resultant force per unit mass, which is denoted in 
the momentum equation by f*. 

(5) Gas is permitted to leave the duct through the walls 
and is denoted in the continuity equation by \|/*, which 
is defined as the mass flow through the walls per unit 
length. 
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It is customary to solve the above set of equations by a numerical 
integration procedure along a set of curves C in the independent variabl 
(x,t)-plane. These curves are usually defined as particle paths and 
characteristic curves. Appendix A presents the derivation of the charac 
teristic equations which are presented here in their final forms. 


S±P 


S In A 

au — gj— - a 


d ln_A + S + S 
dt 5t 


+ " 1 ) a + f 


ta (7-3)/(7-l) e 7(S-So) 


(5) 


5 Q 

^T = ‘ au 


6 S 


^ In A ^ In A , — , ✓ n ^ 

— n a ■ v" — + a + (7 - 1 ) a — - f 

c)x dt 5t v/ ' Dt 


DS 


\|/a 


(7- 3)/f7- 1) 0 7(S-S 0 ) 


( 6 ) 


The entropy condition which must be prescribed for any problem is given 
in a general form as 


DS 

Dt 


F(a,u,S,x, t) 


(7) 


and this completes the system of equations for the three dependent 
variables a, u, and S. The prescribed entropy condition for the cal- 
culation procedures described in the subsequent sections of this report 
is given by 


DS 

Dt 


0 


(7a) 


which characterizes the flow as multi-isentropic. This satisfies the 
condition that the entropy of each gas particle remains constant; how- 
ever, different particles may have different entropies. 


The special symbols for the differential operators are defined as 


+ 



(u + 


\ ^ 
a) ~ 


( 8 ) 


5 



6 


(9) 



+ (u ' *> & 


d s a 

Dt St + U • 


( 10 ) 


All quantities are nondime ns ional in the above equations (see Definition 
of Symbols). Equations (5) through (7) form a system of three linear 
first-order equations for three dependent variables, a, u, and S, that 
will be solved by means of a step-by-step procedure. The parameters P, 
Q, and S vary along curves in the (x,t)-plane that satisfy 


dx _ 
dt “ 

u + a 

for P, 

(U) 

dx _ 
dt " 

u - a 

for Q, 

(12) 

dx 

dt “ 

u 

for S. 

(13) 


The characteristic variables are defined by the relations 

P = a + u (14) 

and 

Q = ~-y a. - u. (15) 


The details of the calculation procedure for solving numerically 
the system of equations (5) through (7) are given in later sections of 
this report. The effects of boundary layer and mass flow through the 
walls of the duct were not considered, and the terms that represent 
those effects have been omitted from the equations. 
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III. the normal shock relations 


The characteristic equations which were derived from the basic dif- 
ferential equations in Appendix A are solved by taking small step-by-step 
increments in time, and the solution for each characteristic proceeds 
along its respective characteristic. These three characteristic curVes 
were referred to as P in direction u + a, Q in direction u - a, and S 
in direction u. Since the characteristic S follows a curve of direc- 
tion u, which is the particle path, two curves o£ this family can never 

cross. Whenever two curves of the same family (either P or Q) meet, a 

discontinuity in the pressure exists at this point. A boundary in the 

flow is thus established, and this boundary is defined as a normal shock 

wave. Two types of shock waves can therefore occur, either a P shock 
(converging of the P characteristics) or a Q shock (converging of the Q 
characteristics). The shock wave path will divide the wave diagram into 
two parts, and on each side certain conditions must be matched. Since 
it is assumed that changes of the flow variables across a shock wave 
take place instantaneously, the steady state relationships between flow 
variables on each side of the shock can be employed. The equations which 
relate the flow variables upstream and downstream of a stationary normal 
shock are generally referred to as the Rankine-Hugoniot equations. The 
derivation of these equations can be found in many references (e.g., 

[12, 13]), and therefore will not be repeated here. Since we are deal- 
ing with shocks that move with respect to the coordinate system, some 
modifications to the stationary shock relations are necessary. The 
equations appearing in this section are presented in their final form 
as they were used in the computational scheme. The shock Mach number Mg 
is defined here as the Mach number of a supersonic flow in which the 
shock would be stationary. 


The relations across a Q shock point are as follows: 



(1) Velocity 

U R • U L „ 2 < 1 ~ M S> 
a L (yfl)Mg 

(2) Speed of Sound 


a R n/ [ 2 +( 7 -l)M|] [2yM|- ( 7 -I) ] 
a L (7+1 )M S 


(16) 


( 17 ) 
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(3) Q Characteristic 


Q R " Q L _ 2 Z^R _ A _ U R ' U L 


a L y ‘ 1 W 


(18) 


(4) P Characteristic 


P R " P L = 2 /^R A U R ' U L 

a T ' 7 - 1 U ' ’ 


(19) 


(5) Entropy 


S R ’ S L 


7(7 - 1) 


In 


{ 


-^-rr M 2 
7+1 S 


Ilk 

7+1 


~> rl + M s 


2±L M 2 J 
2 n S 


( 20 ) 


(6) Pressure 


ll = 27 2 7 - 1 

P L 7 + 1 S 7 + 1 • 


( 21 ) 


(7) Shock Velocity Relative to the Duct 

W S - \ • W (22) 

The subscripts L and R refer to the flow properties on the left- 
and right-hand sides of the shock point, and the subscript S refers to 
the shock. Analogous relations for a P shock point are: 
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(1) Velocity 

U R ' \ 2 d ' M |> 

a R ' (7 + DM S • 

(2) Speed of Sound 

a L _ J[ 2+(7-l)M|][2 7 M|-( 7 -l)] 
^ = (7+1 )Mg 


( 23 ) 


(24) 


(3) Q Characteristic 


Q L " Q R 


7 - 1 \a 


— - 1 + 


U R ” U L 


R 


(25) 


(4) P Characteristic 


p - p 

L R 


2 L 


a R ? ’ 1 W 


— - 1 - 


U R " U L 


(26) 


(5) Entropy 


S L* S R 


~~rr In -j[^rr M? - 


7(7 -1 ) 1L7 + 1 S 7 + 1 J 


r. 


1 + ^ M |-|7- 

^±T^Tj J 

2 S 


(27) 


(6) Pressure 

. ll— 

P R 7 + 1 S 7+1 

(7 ) Shock Velocity Relative to the Duct 


(28) 


( 29 ) 
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The calculation procedure at a shock point involves the matching 
of the characteristic solution with the Rankine-Hugoniot solution at 
the shock point. From equations (18) and (26), it can be seen that 
the relations 


Qo - Q, 


) - (!k IV) 

'0 shock ' S R *\ 


L 'Q shock 


P shock 


(30) 


where for the Rankine-Hugoniot solution, 

2 p/[2+( 7 -l)M|][2 7 M| - (7-1)] ^ 2(1-m|) 

^S R H . = ~ t (7+1 )M g ’ 1 J ' (n l)M g » 

(31) 

and, from the characteristic solution, Asch is calculated from equa- 
tion (30) for the appropriate shock. The subscripts R.H. and CH. refer 
to the Rankine-Hugoniot and characteristic solutions, respectively. 

The condition thatAg R R = Ag CR at a shock point is satisfied by 

expressing equation (31) as 


f (M s ) 



2 fj T2+(7-l)M|][27M| - (7-1)] 
~ \ (7+1 )Mg 


u+ 


2(1 -H|) 

(7+I )Mg 


= 0, 

(32) 


which, with ^s CH and 7 given, is an algebraic equation for Mg. This 

equation is solved in the computer program by the Newton-Raphson method 
for finding the real roots of algebraic and transcendental equations. 
The iteration procedure to satisfy matching the solution of the charac- 
teristic equations with the Rankine-Hugoniot equations at a shock point 
is described in Section VII. 

Some remarks should be made here about how the presence of a shock 
wave (P or Q) affects the flow and possibly the order of calculation as 
pertaining to the wave diagram. Assuming the gas to be flowing from 
left to right, the velocity of a P shock is supersonic relative to the 
velocity of the gas flow to the right; therefore, the P shock overtakes 
the flow of gas ahead of it. The flow conditions which lie to the 
right of the P shock path in the wave diagram are not then influenced 
by either the presence or strength of the P shock wave. The speed of a 
gas which is flowing from left to right will then be increased after 
passage of the P shock wave. 
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For a Q shock, the velocity of the flow is greater than the veloc- 
ity of the shock relative to the duct, and therefore the flow conditions 
which lie to the left of the Q shock path in the wave diagram are inde- 
pendent of the strength or presence of the shock wave. It is seen, then, 
that, in case of a gas flowing from left to right, a Q shock will slow 
down the flow as it passes through the shock. 


IV. THE INITIAL CONDITIONS 




Consider a duct which is closed on the left end and which opens into 
the atmosphere or into an emptying reservoir on the right end divided by 
a diaphragm into two chambers. Let the pressure of the gas on the left- 
hand side of the diaphragm be higher than on the right, and furthermore, 
it is not restricted that the gases on both sides of the diaphragm be the 
same or have the same temperature. When the diaphragm is suddenly 
removed, the gas which was initially on the left side of the diaphragm 
expands and compresses the gas which was initially on the right-hand side, 
forming a P shock. The P shock travels through the gas which was initially 
on the right-hand side. Also created at diaphragm rupture is a contact 
surface. This contact surface is defined as the interface between the 
paths of the gas particles which were initially in contact with each side 
of the diaphragm. The problem at the instant the diaphragm is ruptured 
therefore consists of the simultaneous solution of the characteristic 
equations, the Rankine-Hugoniot equations, and the boundary conditions 
for the contact surface. It is assumed that the diaphragm is instantly 
removed, so the discussion that follows applies for time equals zero at 
the diaphragm location x^. 

Consider the duct shown 
in the illustration which is 
divided into two chambers by 
a diaphragm. Let the sub- 
script 2 refer to the high 
pressure gas on the left- 
hand side of the diaphragm 
and the subscript 1 refer to 
the gas on the right-hand 
side. The prescribed initial 
conditions for the gases in 
the two chambers are the 
following : 


t 
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(1) the pressure of the gases on each side p 2 and p*, 

(2) the ratio of specific heat of the gases / 2 and y x and 

their gas constants R 2 and R x> 

(3) the temperatures of the gases T* and T*, and 

(4) the geometry of the duct in a form such that the term 

(dA/dx)/A can be calculated at any station, x, the loca- 
tion of the diaphragm, xjj), and the total length, L Q . 


It is convenient to choose the initial conditions of the gas on the 
right-hand side of the diaphragm as the reference state. One has then 
the following values for the initial properties of the flow on the right- 
and left-hand sides, respectively. 


a l - -t " 1 


o 

u* 


Ui = 


= 0 


a* 

a Q 


Pi = 


7I^T ai + Ul 



u* 

Up = — 2 = 0 


P 2 “ _ 1 a 2 + U 2 


(33) 


Qi 



Q2 



u 2 




Si - S Q — 0 


s 2 



p 2 -1 

_ a (27 2 )/f7a-l) ' 

2 


All distances are nondimens ionalized by the total length of the duct, 
L q , and the time is nondimens ionalized by 


t 



L 

o 


(34) 
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The wave phenomenon that develops at diaphragm rupture is shown in 
the illustration. The details of calculating the initial solutions are 
given in the following steps: 

(1) Assume that the velocity u 3 (the velocity between the 
contact surface and the P shock) has been given as a result of the pre- 
vious iteration cycle or for the first iteration has been guessed. 

Since the velocity, u 1( and speed of sound, a x , are known, the P shock 
Mach number can be calculated by the quadratic solution of equation (23) 
from the Rankine-Hugoniot relations: 



(35) 


(2) Calculate the flow properties on the left side of the 
P shock by the Rankine-Hugoniot relations; 


a 3 


(7 1 + 1)M« 


• \/[2 + (jl - 1)M|][27 1 M| - (y x - 1)] 


(36) 


= 


7i(7i - 1) 


In 


t 


iZi. M 2 . ZjfcJL. 1 

7i + 1 S 7x + 1 


rl + ^2 Mg-i^ 1 ^. 


21 ± 1 M 2 

2 n s 


} 


(37) 


£3 = 2 ^1 M 2 _ 7 _± ' 1 

Pi 7 1 + 1 S 7i + 1 


(38) 


(3) One of the boundary conditions across the contact surface 
is satisfied by setting p 4 = p 3 , and since the expansion takes place 
isentropically, the speed of sound on the left side of the contact sur- 
face is given by 


a 4 ” a 2(p 4 / P 2 ) 


( 72 ‘ 1 )/ 2 72 > 


(39) 
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(4) The velocity u 4 on the left side of the contact surface 
is given by the requirement that the characteristic variable P remains 
constant from the head-wave to the left side of the contact surface: 


P 4 = P 2 . 

From the definition of P, the velocity u 4 can be calculated by 


( 40 ) 


u 4 



(41) 


(5) Satisfy the second boundary condition across the contact 
surface by setting u 3 = u 4 which yields a revised value for u 3 . Steps 
(1) through (5) are repeated until the difference between the calculated 
values in u 3 for subsequent iteration cycles is less than a prescribed 
tolerance. 

At the initial time t = 0, a centered expansion fan originates at 
the diaphragm location which is bounded on the left by the head-wave 
(Q 2 characteristic) and on the right by the tail-wave (Q 4 characteristic). 
Additional Q waves are introduced between the head- and tail-waves to 
keep the net of characteristics sufficiently fine; e.g., a total of n 
initial Q waves is employed in the computer program. 


V. CHARACTERISTIC SOLUTION 


In the numerical solution of the set of characteristic equations 
given in Section II, several procedures could be employed. The under- 
lying process for solving these equations is, however, the same regard- 
less of the particular procedure and is the integration of the 
characteristic equations along their respective paths in the x,t-plane. 
This integration is carried out by a step-by-step procedure along short 
straight line segments, the slopes of the respective characteristic 
variables being approximated by the average of their values at the end 
of the short line segments. The characteristic curves for the case at 
hand are designated as P, whose path in the x,t plane is in direction 
u + a, Q in direction u - a, and S in direction u. After exploring 
several different procedures, it was decided to use a constant time 
interval and to follow the path of one of the family of characteristics 
(Q), and to use an interpolation technique to follow the path of the 
other family (P) between subsequent time intervals. The order of cal- 
culation as performed by the computer program is from left to right 
along a line t = constant. 
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A. Calculation Procedure for a Regular Point 


Assume that a distribution of n points has previously been computed 
or given as initial conditions on the line t = t Q and that all properties 
of the flow are known at these points . 

Let it also be assumed that the 
solution has been completed for 
j points on the line t x - t 0 + At 
(see illustration). Let the sub- 
script A refer to properties at 
the j + 1 point on the line t x 
which is to be calculated, B 
refer to the point on the t 0 line 
through which the P character- 
istic that intersects points A 
and B passes, C refer to the 
point on the t 0 line through 
which the entropy path that 
intersects points A and C passes, 

and D is the origin of the Q characteristic from the t Q line which 
passes through point A. The calculation procedure consists of an 
iteration process for the location and properties at point A. 

The straight line segment for the P wave that intersects points A 
and B in the time interval t-)_ - t 0 = At must satisfy the relation 



X A " X B (U + a) A + (U + a) B 
At 2 


(42) 


which can be rewritten as 


At , . , At , . 

X A ' T (u + a) A *B + T (u + a V 


(43) 


The left side of the above equation relates the position and flow pro- 
perties u A and a A at point A on line t x to the corresponding properties 
at point B on the line t Q . Keeping in mind that the flow properties 
are a known function of x on the line t as a result of the previous 
stage of calculation or as given initial conditions, then for given 
quantities of x^, u A , and a^, we can find the point xg and the flow 
properties at Xg by interpolation if a table of 

z x (x) = X + (u + a) (44) 


has been prepared for the distribution of points on the t Q line. 
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Following the same approach for the entropy path used for the P 
wave, the straight line segment for the entropy path in the time interval 
must satisfy the relation 


u A + u 
A c 


At 


(45) 


which can be written as 


At 

X A ' 2 U A 



(46) 


A table of values for 


z 2 (x) ■ x + y u (47) 

is also assumed to have been prepared and stored at the completion of 
the calculations for the t 0 line. Since the procedure follows each Q 
characteristic from one time interval to the next, an analogous table 
for the Q wave is not necessary. 

The step-by-step procedure for calculating the flow properties and 
position of point A is as follows: 

(1) Assuming that u^ and a^ have been given from a previous 
step, or for the first iteration have been guessed, the location of x^ 
is given by intersection of the line segment for the Q-wave path that 
passes through Xp with the constant time line, t = t-^ 

x A = x Q + “■ [(u - a) A + (u - a) D ]. (48) 


(2) Calculate the left-hand sides of equations (43) and (46): 


‘1 ■ *A ' f <» + *>A 

(49) 

At 

C2 = X A • T V 
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whereby the points xg and xq can be found by an inverse interpolation 
of the prepared tables of z x (x) and z 2 (x). The flow properties u, a, 
and S are then given for these two points by interpolation. Values for 
the characteristic variables P and Q are calculated by equations (14) 
and (15) for these two points. 


by 


(3) The characteristic variables P A and Q A are then calculated 


P. = 


P B + 


/ 1 dAv , , 1 dAx 

r(- au 7 — ) 4* (-au 7 — ) 

1 A dx B A dx A 


At + ^ < S A ' V 


(50) 


= + 


, dA\ , , 1 dA\ , 

( ' aU A d^D + ( ' a " A a ' + 3 


At + <S A - S D , 


where the 


1 dA 

geometrical term 7 — is assumed to 
A dx 


be a known function of x. 


(4) From the definition of the characteristic variables, equa- 
tions (14) and (15), new values of the velocity u^ and speed of sound a^ 
can be calculated. 


“A - 2 < P A - V 


“a " 2 "*" 1 (P A + V 


(52) 


(53) 


(5) Steps (1) through (4) are repeated until the differences 
of u^ and a^ for subsequent iteration cycles are less than a prescribed 
tolerance. 

After the convergence is met, the final values of the flow properties 
from the last iteration cycle are stored for use in calculating the solu- 
tions at the next time interval. All interpolations that were carried 
out in the above steps are linear. 
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B. Calculation Procedure for Points Near Boundaries 


The procedure for calculating a regular point given under Section V(A) 
involved projecting a Q wave from time t Q to time t-L along its respective 
characteristic, and by an interpolation process finding the point at the 
previous time interval where the P characteristic passes that intersects 
with the Q wave at the time t x . It is obvious that this standard procedure 
cannot be applied for calculating points very near to the right side of a 
boundary or discontinuity since the Q wave and entropy path can cross the 
boundary in the time interval. For example, the illustration shows the 
intersection of the P and Q characteristics along with the entropy path 
at point A. Since the P characteristic intersects the Q shock at point 
B, the standard procedure given in Section V(A) is not applicable. 



An alternate procedure was used for calculating points on the right 
side of a boundary .where the P characteristic or entropy path crosses 
the boundary in the time interval under consideration. Consider the 
path of the boundary between times t Q and ^ along with the flow pro- 
perties on the right-hand side of the shock at times t Q and t x to be 
known. The equations that must be solved to determine the location and 
flow properties at point A are: 



( u + a) A + (u + a) B 
2 


along the P wave, 


( 54 ) 
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along the Q wave, 


(55) 


X A - *D (u ~ a> A + (u ' a) D 
At 2 


X A : *C U A + U c 

t, - t 2 


S A = S C> 


along the entropy path, 


(56) 

(57) 


P = P + 
A B 


/ 1 dA\ . , 1 dAs 

(-au — -r— ) + (-au — -t— ) -| 
A dx'A A dx y B 1 


(t, - v 


a + a 

H — h: ® . 3 ) 

2 ^A ®B' ’ 


(58) 


^A = Q D = 


, 1 dAv , , 1 dA, 

aU A dx^A ^ aU A dx^D^ 


a + a 

At +-S- E <S a - V' 


(59) 


The equation for any point B on the path of the boundary between 
times t Q and t x is given by 


x 


R, 1 ' *B 



(ti - 


V* 


(60) 


where the subscripts R,0 and R,1 refer to properties on the right-hand 
side of the boundary at times t 0 and t 1# Combining equations (54) and 
(60) results in 


t-, - t = 


x - 

A R, 1 


(n + a) A+ (u + a) B x R l - x r>q 


At 


(61) 


which relates the intersection of the P wave with the boundary. Replac- 
ing the subscript B in equation (60) by the subscript C and combining 
this equation with equation (56) results in 
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I 


(62) 


t;L " t C u, + u 


x A - x„ , 
A R.l 


C _ X R,1 " X R, 0 


At 


which relates the intersection of the entropy path with the boundary. 

The step-by-step procedure for solving the above set of equations 
for the position and flow properties at point A is as follows: 

(1) Assume that u^ and a^ are known from the previous itera- 
tion step or for the first iteration have been guessed. Equation (55) 
gives the relation for x^ 


X A = *D 



[(u - a) A + (u 


a) D l* 


(63) 


(2) The point of intersection of the P wave and the boundary 
is found by solving equation (61) by the iteration process, i.e.. 



(a) 

guess t B , 


(b) 

obtain Ug and ag by a linear interpolation of 
these known properties on the right side of 
the boundaries at times t Q and t x , and 


(c) 

calculate an improved value of tg from equa- 
tion (61). 

Steps (b) and (c) are repeated until changes in tg satisfy a prescribed 
tolerance. The axial location xg is then calculated by equation (54), 
and Sg is given by interpolation as was done for ug and ag in step (b). 


(3) The point of intersection of the entropy path and the 
boundary is given by solving equation (62) by the same method as 
described above for equation (61). The entropy S c is then given by 
interpolation of the known values of S at x-^ R and Xq r . The entropy 
at point A is then given by (57). If t c § t^, the entropy is handled 
in the same way as described in Section V, paragraph A, for a regular 
point. 


(4) Values for the characteristic variables at point A are 
given by equations (58) and (59) from which improved values of u A and a^ 
can be calculated. 
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(64) 


a A = L T 1 (P A + V 


U A 2 (P A " V' 


(65) 


Steps (1) through (4) are repeated until u^ and a^ converge to a 
final value within a prescribed tolerance. 


VI. CALCULATION PROCEDURE FOR A CONTACT SURFACE 


A contact surface is defined as the interface between the paths of 
two different gases, or the same gas at different entropy levels, in 
which no flux of matter passes. In the x,t-plane the contact surface 
is a line of discontinuity in which the speed of sound and the entropy 
are discontinuous. By nature of its definition, the flow velocity and 
pressure through the contact surface at a point in the x,t-plane are 
cons tant . 

Before entering into a discussion of the iteration procedure for 
calculating the flow properties on each side of the contact surface, 
some relations which will be required are derived. Let the subscripts 
L and R refer to flow conditions on the left and right sides of the 
contact surface, respectively. The characteristics P on the left and 
Q on the right are, by definition, 


P 


L 



a L + 


U T 


( 66 ) 


Q 


R 


• V 


(67) 


By applying the boundary condition u L = u R , equations (66) and (67) 
combine to give 


P T + Qr 



a L + 



V 


( 68 ) 
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The pressures on each side of the contact surface are expressed as 


2 7 l / (7 L “ 1 ) -7l s L 
P t = a T e 


( 69 ) 


27 r /(7 r “ 1) -7 r S 
P = a„ e 

* D 


(70) 


which includes the condition that the reference entropy level is zero. 
Applying the boundary condition p-^ = p R> equations (69) and (70) combine 
to give 


a 


R 





( V 1)/27 r 


(71) 


When equations (68) and (71) are combined, the satisfaction of both 
boundary conditions across the contact surface is included. We then 
obtain 


7r -1 


P + Q = — 
L V R 7 t 


T a + — 

1 L 7 r 


- 1 a L 


[7 l ( 7 r - 1 )^r(7l- 1 )] 2 7 : 


( 7rV 7 l s l> 


(72) 


which, for P L , Q R , S R , S L , y L , and y R given, is a transcendental equation. 

This equation must^be solved numerically for a^. For 7 l = 7r = 7? can 
be solved directly to give 




21 f 1 < P L + 


X 


- 1 


1 + e 


(S R - V 


(73) 
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It is assumed that the solu- 
tion is given for a distribution 
of n points on the t D line and 
that a contact surface passes 
through this line. Also assumed 
is that j regular points have 
been calculated on the t-L line. 
The procedure for calculating the 
contact surface point on the t x 
line is given in the following 
steps . 


n 

Q Waves 
/ 



d 



/ Contact 

/ ^ Surface 

/ 

^ 

X 

0, L * 0 ,R 

> 


(1) The values for the entropy on each side of the contact 
surface are given since the entropy level for each side remains con- 
stant throughout the trajectory of the contact surface in the x,t-plane; 
therefore, set 


S 1,L 4 S 0,L 
S 1,R " S 0,R’ 

where the first subscript 0 or 1 refers to the t Q or t x line, respectively. 

(2) Assuming that u^ L and a^ L have been given from the pre- 
vious iteration cycle, or for* the first iteration cycle, have been 
guessed, calculate the x location of the contact surface point on the 
t ± line: 


X 1,L X 1,R 




a 0,I? 


(74) 


(3) Satisfy the boundary condition of no velocity change across 
the contact surface at a point by setting u^ ^ = u^l and satisfying the 
analogous boundary condition for pressure by calculating the speed of 
sound on the right-hand side from equation (71). 

(4) Find the points on the t Q line through which the P and Q 
waves that intersect at x^ L = x 1 R pass. This is accomplished by an 

interpolation process that’ locates the x stations on the t Q line that 
satisfy 
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( 75 ) 


X 1.L ' X P (u + a >l,L + (U + a) P 
At 2 

for the P wave, and 


x 


LR 


At 



(u - a) ltR + (u - a ) Q 
2 


(76) 


for the Q wave. The subscripts P and Q refer to properties at the points 
xp and xq on the t Q line which satisfy the above relations. 

(5) Calculate values for l and R by way of the character- 
istics equations 


P = P + 
1,L P 


, 1 dA. , . 1 dAv 

+ ( ' a " A teh.L 


At + 


P + a l,L 


and 


Vr = % + 


/ 1 dA \ , / 1 dA \ 

( ' a “ A d7>0 + ( ~ aU A SYr 




(S. t' S o> 

1,L P 
(77) 


(S 1,R-V- 


(78) 

(6) Calculate an improved value for a^ ^ by solving equation 
(72), or by solving (73) for _ y R , and an improved value for u-^ L by 


U 1,L P 1,L “ 7 l - 1 a l ,L‘ 


(79) 


Steps (2) through (6) are repeated until both u-^ ^ and a-^ ^ converge to 
constant values within a prescribed tolerance. 
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VII* CALCULATION PROCEDURE FOR A SHOCK POINT 


As suggested by Hartree [6], it was found that the treatment of 
boundaries and discontinuities was simpler when using a constant time 
interval rather than following a grid of characteristics. This method, 
explained in detail in Section V for a general point, follows one of 
the family of characteristics (Q waves) from time t G to time = t Q +At 
and interpolates the x station at time t 0 through which the other 
characteristics (P waves) pass. The order of the calculations proceeds 
from left to right in the wave diagram for each line t * constant. 
Because of the order and method by which the calculations are carried 
out, the detailed procedure used to calculate a Q shock point is some- 
what different from that for a P shock point. The details of the pro- 
cedure for each type of shock point will therefore be presented 
separately. 


A. Q Shock Point 

It is assumed that the solution has been completed for a distribu- 
tion of n points on the line t = t Q and that a Q shock passes through 
this line. Since changes of the flow variables across the shock are 
assumed to take place instantaneously and since the shock thickness is 
assumed to be zero, the properties of the flow at the shock point are 
double-valued. To explain the procedure for calculating a Q shock 
point, a double subscript of notation is used. Let the first sub- 
scripts 0 and 1 refer to the shock at times t 0 and t Q + At and the 
second subscripts L and R refer to the left- and right-hand sides of 
the shock, respectively. It is also assumed that j regular points 
which follow Q characteristics have been calculated on the t Q + At 
line where the j-th point is the 
point generated from the point on 
the immediate left side of the 
shock wave on the t Q line. If 
we keep in mind that Q charac- 
teristics overtake a Q shock 
from both sides, the illustra- 
tion shows Q waves that have 
crossed the shock path from the 
left side. After the iteration 
for the shock point has been com- 
pleted, those points that lie on 
the right-hand side of the shock 
at time t x are discarded. These 
points were retained only for 
interpolation of properties on 
the left side of the shock dur- 
ing the iteration procedure for 
the shock point. 


t 



x 

25 



The calculation procedure for the Q shock point on the t x line is 
now given in the order that the computer program actually solves the 
problem: 


(1) Guess the shock velocity w^ g at time t x = t Q + At. 
Calculate the average shock velocity between times t Q and t^: 


w = 


"o.s +W 1.S 


(2) Calculate the position of the shock point on the t ± line: 


= x. 


1 , L 1 , R 


■ *0,L + At 


w. 


(80) 


(3) Interpolate for the properties u^ L , a^ L and L on 
the left side of the shock point. 

(4) Calculate the shock Mach number: 


M 


1,S 


U 1,L ~ V 1,S 


1,L 


(81) 


(5) Calculate the properties on the right side of the Q shock 
point which are given by the Rankine-Hugoniot equations presented in 
Section III. 

(6) The point Xq on the t Q line through which the Q character- 
istic passes that also passes through the point x-. R is found by satis- 
fying 


-i _ — x 

1 > R Q. = 

At 


(U " a) l,R + 
2 



(82) 


This equation is solved by the 
method of iteration for x^. 
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(7) Calculate a new value for the Q characteristic on the 
right side of the shock point by 


Vr q q + 


/ 1 dA \ . f 1 d A \ 

^~ aU A dx^Q ^~ SU A dx^l ,R 


^ + -¥- 2 (S 1 ,r-S q )- (83) 


(8) The following relationship across the shock is now 
calculated: 


^ = 




(84) 


This relationship yields a new value for the shock Mach number g 
through the Rankine-Hugoniot equation (32). * 


(9) A new shock velocity w^g is then given by 


W 1,S U 1,L ” a l,L M 1,S' 


(85) 


A comparison of the guessed shock velocity with the calculated shock 
velocity determines whether a prescribed tolerance in their difference 
has been met. Steps (1) through (9) are repeated by replacing the 
guessed value of w^s by the calculated value until the tolerance is met. 

All points, previously calculated on the line t-,_, which lie on the 
right-hand side of the shock point are now dropped. Calculation of points 
on the right of a boundary is discussed in Section V, Paragraph B, of 
this report. 
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B. P Shock Point 


The calculation procedure for a P shock point is somewhat different 
from that of a Q shock point, but the underlying process of matching the 
Rankine-Hugoniot solution with the characteristic solution at the shock 


point is the same. Again it is 
assumed that the calculation pro- 
ceeds from left to right and that 
j points have been calculated on 
the t x line as shown in the 
illustration. The detailed pro- 
cedure for calculating the P 
shock point is given in the 
following steps. 
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(1) Guess the shock velocity 5 and calculate the average 
shock velocity between times t Q and t ± : 

- ”o.s + "l.S 

w = — 2 — 2 • 

(2) Calculate the position of the shock point on the t ± line: 


k l,L 


= X 1 , R = X 0,L + At ' W - 


( 86 ) 


(3) Find the locations on the t Q line through which the P and 


Q characteristics pass that intersect at the point R . 
plished by approximately the 
same process as that for cal- 
culating a general point, and 
therefore a step-by-step explana- 
tion of the procedure will not be 
repeated here. The difference is 
that a particular Q characteristic 
is not followed; therefore, inter- 
polations at time t Q are neces- 
sary to find both the P and Q 
characteristics that pass through 
the point xx r. This step yields 
the flow variables on the right- 
hand side of the P shock. 


This is accom- 



x 


(4) By knowing the flow variables on the right side of the 
shock point from the above step, the P shock Mach number is calculated 
by 


w 


M 


1 iJL 


jLR 


i,s 


1,R 


(87) 


(5) The Rankine-Hugoniot relations (Section III) give the flow 
variables on the left-hand side of the P shock. 


( 6 ) The point xp on the t Q line through which the P character- 
istic passes (see illustration) and that also passes through the point 
xx l f° unc l by satisfying 


x i.l : x p 

At 


(u + a)p + (u + a) L L 

2 


This equation is solved by the method of iteration for Xp. 
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( 88 ) 



(7) Calculate a new value for the P characteristic on the left- 
hand side of the P shock by 


* 


1 >L 


P + 
P 


(*au 7 ^ ) + 

A dx P 




At + 


l.L +a P 


(S 1,L ’ V 


(89) 


(8) Calculation of the relationship 


= -LL Lil (90) 

a l,R 

across the shock yields a new value for M]^ 5 by way of the Rankine- 
Hugoniot equation (32). 

(9) A new shock velocity w-^ g is calculated by 


w. ^ = u -« + a -. „ M- 

1,S 1,R 1 ,R 1,S 


(91) 


Comparison of the guessed value with the calculated value of the shock 
velocity determines whether a prescribed tolerance has been met for their 
difference. Steps (1) through (9) are repeated by replacing the guessed 
value with the calculated value until this tolerance is met. 


VIII. INTERACTION OF DISCONTINUITIES 


Thus far in this report, we have dealt only with cases where charac- 
teristics cross or intersect; however, it can easily be seen that there 
are many other possibilities of boundaries and discontinuities intersect- 
ing with one another, such as shocks of like families intersecting, shocks 
of unlike families intersecting, and shock contact surface intersections. 
Although the general methods for handling these interactions are similar, 
each case requires somewhat different calculation procedures. Because of 
the large number of possibilities that can occur, only the intersection of 
the discontinuities that were encountered in the problem at hand will be 
discussed. The reader should be able to easily modify the discussed pro- 
cedure for a particular interaction not treated herein. 
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A. Merging of Two Q Shocks 


When one Q shock overtakes another, the wave diagram in the vicinity 
of the point of intersection of the two Q shocks appears as shown in the 
illustration. A Q shock stronger 
than either of the merging shocks 
results. Also produced is a con- 
tact surface which is the par- 
ticle path that passes through 
the intersection point and 
separates the gases that have 
been compressed by the two 
shocks. For 7 ^ 5/3 (see ref- 
erence 10), a reflected wave 
which was created by the inter- 
actions was found to be quite 
weak for all of the shock- 
shock interactions encountered 
in this problem. 



The interaction of two Q shocks was solved by an approximate trial- 
and-error procedure. Assume that two Q shocks have been detected to 
cross between times t Q and t x . Of course, there must be some logical 
check to detect this crossing in the computer program, but, as there 
are numerous ways this could be handled, this will not be discussed 
here. Let tj refer to the point of intersection of the two Q 
shocks, the subscripts 1 through 5 refer to flow properties of the 
regions in the vicinity of the interaction as shown in the illustra- 
tion, and A, B, and C refer to the three Q shocks. The iteration pro- 
cedure is as follows: 


(1) Consider that the velocity u 5 has been given from the 
previous iteration step or for the first iteration cycle has been 
guessed. It is assumed that the flow properties of Region 1 are given 
by the values of the flow which had been calculated on the left side of 
shock A at time t Q , and that the property values of region 3 are those 
which were calculated for the right-hand side of shock B at time t Q . 
These assumptions are justified because of the smallness of the time 
step At. 


(2) Properties in region 5 can then be calculated from the 
Rankine-Hugoniot relations across a Q shock: 


Ai 


C 


- u 


(92) 
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2 


. z 


+ 1 


+ 


M c = 


2 + 1 


Au, 


+ 4 


(93) 


a 5 = a ; 


\/[2 + (7 - 1)M 2 ] [2yM 2 - (7 - 1) 

(7 + 1)M_ 


(94) 


S= = S i + 


In I M 2 - — 

7(7 - 1) in l]_7+l n c 7+1 _ 


rl +^ 


M" 


M 2 J 
2 C 


7) 


(95) 


(3) Since the reflected wave was a centered expansion fan for 
the range of 7 possible for this problem, the entropy across the expan- 
sion is constant; i.e., S 4 = S 3 . The speed of sound for region 4 is 
then given by satisfying the boundary condition of equal pressure across 
the contact surface: 


a 4 - a 5 e 


— (S4-S5) 


(96) 


(4) The reflected wave, which is always an expansion wave for 
7 ^ 5 / 3 , separates region 3 from region 4, and the value of the charac- 
teristic variable Q across the reflected wave is constant. Therefore, 
by setting Q 4 = Q 3 , the velocity in region 4 is given by 

u 4 = _ i a 4 - Q 4 . (97) 

(5) An improved value for u 5 is given by satisfying the second 
boundary condition across the contact surface: 


u 5 = u 4 . (98) 

Steps (2) through (5) are repeated until the difference in calculated 
values of a 5 for subsequent iteration cycles is less than a prescribed 
tolerance . 
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The shock velocity of shock C is then given by 


w c = ui + a 1 M c> 


( 99 ) 


the location of the Q shock point on the line t x is 


x c (t i) = Xj + w (; (t 1 - t ]; ), 


( 100 ) 


and the location of the contact surface point on the line t x is 


x c.s. (tl) = x i + U5(tl ■ t i ) ' 


( 101 ) 


The reflected wave was found to be very weak for all the interactions 
of two Q shocks that occurred in the problem. Points on the right-hand 
side of the contact surface on line t x were calculated in the same manner 
as the points near boundaries (see Section V, Paragraph B). This, in 
effect, spreads the centered expansion, which was weak over the time 
interval At rather than at a point. This approximate way of handling 
the reflected wave should be quite good since the time step is small. 


B. Crossing of a Q Shock by a Contact Surface 

The following illustrations show the two types of interactions that 
can occur when a contact surface intersects with a Q shock. 
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If the pressure ratio across the transmitted Q shock p 5 /p 2 is less 
than the pressure ratio across the incident Q shock P 3 /p 1? then the 
reflected wave is a rarefaction fan. However, if p 5 /p 2 > p 3 /pi, the 
reflected wave is a P shock. In the cases calculated for the problem 
at hand in which a contact surface crossing a Q shock occurred, the 
reflected wave was always a weak rarefaction fan. The computer program 
was coded for the case of only the reflected wave being a rarefaction 
fan, but maintaining the capability of detecting if the reflected wave 
is a P shock and printing out a message describing this condition and 
the strength of the P shock. The rarefaction fan is handled in the same 
manner as that described in the previous section when two Q shocks 
interact. 

The interaction of the contact surface and a Q shock was solved by 
an approximate trial-and-error procedure. Let xj, tj refer to the inter- 
section point, the subscripts 1 through 5 refer to the regions indicated 
in the illustrations, and A and B refer to the incident and transmitted 
Q shocks, respectively. Flow properties in regions 1, 2, and 3 are 
assumed known. Calculation steps for solving the interaction process 
are as follows: 

(1) Consider that the velocity in region 5 (u 5 ) has been given 
from the previous iteration cycle or for the first cycle has been guessed. 

(2) The remaining properties in region 5 can then be calculated 
from the Rankine-Hugoniot relations across a Q shock: 


M. 


B 



Au b 



2 


(103) 


a 5 - a 2 


n/[2 + (7 - 1)M|] - (7 - 1)] 

(7 + 1 )^ 


(104) 


s 5 - s 2 + 


7(7 - i) 


In 


-P 2_ „ 2 . izl 


|]_ 7+1 B 7 + 1 J I 2± i. m 2 

2 B 


rl +z r M f- 


n 

I 


(105) 


33 


(3) Since it has been assumed that the reflected wave is a 
rarefaction fan, there is no change in entropy across the reflected wave; 
i.e., S 4 = S 3 . The speed of sound is then given by satisfying the bound- 
ary condition of equal pressure through the contact surface. 


Zli 

2 

a 4 ” a 5 e 


(S4-S5) 


(106) 


(4) An additional condition that must be satisfied is the con- 
stancy of the characteristic variable Q across the expansion fan 


Q4 ~ Q3 


which by definition yields the relationship for the velocity in region 4. 


U4 = — 1 a 4 - Q4 • 


(107) 


(5) The second boundary condition across the contact surface is 
now satisfied by setting 


u 5 = u 4 


which then provides the improved value for u 5 in the next iteration cycle 

Steps (2) through (5) are repeated until the prescribed tolerance 
for the difference in successive values of u 5 is met. 

The location of the shock and contact surface points, and the cal- 
culation of points on the right side of the contact surface by the charac 
teristic solution at time t 1? follow the same procedure described in the 
previous section when two Q shocks interact. 


IX. EXIT CONDITIONS 


After diaphragm rupture, the volume of gas that has been compressed 
on the left-hand side of the diaphragm flows through the facility and 
empties into a reservoir which is connected to the end of the duct. 
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Since the mathematical solution is intended to calculate only the early 
development of the flow, it is not necessary to include the conditions 
for inflow at the end of the duct. The possibility of inflow at the 
exit could exist only at a very late time because of a buildup of pres- 
sure in the emptying reservoir, or at a very late stage of the blowdown 
if emptying into the atmosphere. 

For the problem at hand, gas is always assumed to be leaving the 
duct at the exit station L Q and is therefore treated as an outflow 
problem. The gas in the external region is assumed to be at rest. The 
calculation procedure for satisfying the appropriate boundary conditions 
for subsonic, sonic, or supersonic flow at the exit is now given. 


A. Subsonic Flow at Exit 


For subsonic flow it is assumed that the pressures on each side 
of the exit plane follow the steady state boundary condition of equal 
pressure. Since we have chosen the conditions which were on the right 
side of the diaphragm before rupture as the reference state, the pres- 
sure at the exit is given by p EL = p Er = 1 where the subscripts E E and 
Er refer to the left- and right-hand sides of the exit plane, respec- 
tively. The iteration procedure to calculate the flow properties at 
the exit is as follows: 

(1) The velocity u^ is given by the previous iteration cycle, 
or for the first cycle is guessed. 

(2) The entropy S E at the exit is given by the conventional 
procedure as for a regular point (Section V, Paragraph A). 

(3) Applying the boundary condition p = 1 at the exit results 
in the following equation for the speed of sound, we obtain 


a E 


L 



which includes the conditions that the properties on the right side of 
the exit plane are the reference state and the reference entropy is zero. 

(4) The value for the P characteristic P E at the exit is 
found in the same manner as for a regular point. This value then fur- 
nishes a relation for the velocity: 
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Steps (2) through (4) are repeated until the change in u E for sub- 
sequent iteration cycles is less than a prescribed tolerance. L 


B. Sonic Flow at Exit 


The boundary condition to be satisfied for sonic flow at the exit 
is u E , = a E ^. The tr ial-and-error procedure to calculate the sonic flow 
conditions at the exit is as follows: 

(1) Assume that the velocity ujr^ has been given from the 
previous iteration cycle or for the first cycle has been guessed. 

(2) Satisfy the boundary condition by setting a E ^ ■ u E ^. 

(3) Values for the entropy Sg L and characteristic P Ee are found 
by the standard procedure (Section V, Paragraph A). 

(4) A new value for u E ^ is given by 


= 21 


+ 1 % 


Steps (2) through (4) are repeated until the change in u E ^ for sub- 
sequent iteration cycles is less than a prescribed tolerance. 


C. Supersonic Flow at Exit 

For supersonic outflow, the exit conditions can be calculated by 
the standard procedure since both P and Q characteristics reach the exit. 


X. SOME PERTINENT DETAILS ABOUT THE CALCULATION PROCEDURE 


A. Geometrical Aspects 

Effects on the flow properties due to the change in cross section 
are felt through the term 

I dA 

A dx 


which appears in both characteristic equations. Information concerning 
the geometry of the duct must be furnished to the computer program in a 
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form such that the term 


I dA 

A dx 


can be evaluated at any point x. The actual duct shape is approximated 
by a series of straight line segments for the inside diameter in the 
axial direction. To avoid discontinuities in the first derivatives of 
the area at the junction of the line segments, a parabola which con- 
nected the straight lines at a distance of e on each side of the 
junction was used as shown in the illustration. 


The computer program requires M input 
cards regarding the geometry for which 
each card contains an axial station x^(i) 
along with the corresponding diameter 
d|(i), and curve fit parameter e|(i), where 
i = 1 ,M. The left end of the supply tube 
x* = 0 is referred to by i = 1, and the 
right end of the facility where x* = L 0 
is referred to by i = M. The remaining 
indices i = 2,M-1 refer to the junction 
of the straight line segments. The program 
automatically nondime ns ionalizes all dis- 
tances by L 0 . 



B. Grid Control 


To keep the characteristic grid at a sufficient density, it was 

necessary to establish some criteria to either introduce or take out 

points at each time step. It is seen, for instance, that in regions 
where there is sonic flow, additional points must be introduced, since 
in subsonic flow, the Q characteristics travel to the left, and in 
supersonic flow, they travel to the right in the wave diagram (see 
Figure 2). 

Whether or not points should be introduced was determined by pre- 
scribing an upper limit in the change of Q between two consecutive 
points and an upper limit in the allowable distance between them. Let 
these limits be referred to as an< ^ ^max anc * consider that the 

solution has been completed up to a point j on a line t = t x . The 

following quantities are then calculated: 

= Q(j) - Q(j-l) 

Ax = x(j) - x(j-l) . 
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If either A} or Ax is larger than the prescribed upper limits for these 
quantities, any number from one to three additional points or Q waves are 
introduced between x(j) and x(j-l). These additional points are calcu- 
lated by the method of characteristics in a special subroutine called 
POINTADD of the computer program. The maximum allowable difference in 
Q between two adjacent points £Qntax is a required input (TOLQ) of the 
computer program. Values of TOLQ =.05 for the subsonic cases and 
TOLQ = .2 for the supersonic cases were found to give good results. 

The upper limit for the distance between any- two points AXmax is auto- 
matically determined by the computer program. The value of Axmax is not 
constant throughout the duct since, in sections where there are area 
changes, a much denser distribution of points is necessary than in 
straight sections where there are no area changes. The program also 
takes out points automatically when the grid becomes too dense. 

By running repeated cases for several time increments, it was found 
that At = .001 for subsonic and At = .002 for supersonic cases gave 
satisfactory results. What is meant by subsonic and supersonic cases 
as referred to above pertains to the state of the flow in the test sec- 
tion after the flow is established. The distance between the points in 
the axial direction was found to be much more critical than the time 
interval in the sections where there were large area changes. 


C. Restart Procedure 

Since each case took a fairly large amount of time on the CDC 3200 
computer, a restart procedure was incorporated in the program. Each 
time on the computer, a prescribed number of time intervals are computed 
and the distribution of points and their pertinent properties for the 
last calculated time are punched on computer cards, which are then used 
as input for the next run. The preparation of input for the program 
is explained in detail in Appendix B. 


XI. DISCUSSION OF THE RESULTS 


The calculation procedures described in this report for the one- 
dimensional flow through the pressure- tube tunnel were programmed for 
the CDC 3200 computer. Calculations were performed for the anticipated 
range of Mach numbers for the facility. Figure 1A presents the aero- 
dynamic boundaries for the full-scaled facility with the Mach 2 nozzle 
installed. The interchangeable nozzle sections for the supersonic cases 
were designed from the results of a two-dimensional method of character- 
istics satisfying the condition of uniform flow at the nozzle exit plane. 
The spool length satisfies the condition that the length of the nozzle 
and spool equaled 138 inches for each case. The distance from the 
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diaphragm to the nozzle throat was approximately 29 percent longer for 
Mtest = 3.5 and 5.0 than the corresponding length for M<eest = 1*4, 1.7, 
and 2.04. Figure IB presents the aerodynamic boundaries of the tunnel 
with the sonic nozzle installed. This configuration was used for the 
sonic and subsonic cases that were calculated. The subsonic Mach 
numbers in the test section were attained by choking the flow downstream 
of the test section by the use of choking flaps as shown in Figure IB. 
Also shown in Figures 1A and IB are imaginary or effective boundaries 
which were assumed to approximately account for blockage by the model 
support and diaphragm cutter mechanisms. 

The following initial and reference conditions were assumed for all 
of the cases calculated: 

(1) The gas on each side of the diaphragm is air whose ratio 
of specific heats is y = 1.4. 

(2) The temperature of the air on each side of the diaphragm 
is T* = T* = 295.57 °K. 

(3) The properties of the air on the right side of the 
diaphragm before rupture are: 

a'^ = 1130 ft/sec 

p* = 14.7 lb/in 2 = 2116.8 lb/ft 2 
S x = 0. 

(4) The reference conditions are the state of the air on the 
right side of the diaphragm before rupture. 

(5) The reference length is the total length of the facility: 

L q = 325 ft. 

Figure 2 presents the wave diagram for the case in which the Mach 
number in the test section during the period of steady flow (referred 
to as M^est) was two - Briefly recapping what occurs during the early 
development of the flow, we can see the centered expansion fan, which 
is bounded on the left by the head-wave and on the right by the tail- 
wave, that originates at the diaphragm location at the instant the 
diaphragm ruptures. The expansion waves that make up this centered 
expansion fan accelerates the air which was initially on the left side 
of the diaphragm so that it starts to flow out through the facility. 

Also shown is the P shock that was formed when the diaphragm was 
ruptured. The P shock accelerates the still air through which it passes 
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as it travels out of the duct. The path of the contact surface, i.e., 
the interface that divides the two volumes of air of different entropy 
that were initially in contact with each side of the diaphragm before 
rupture, is shown as it travels out of the duct. Between the tail-wave 
of the initial centered expansion fan and the contact surface, two Q 
shocks are formed which eventually interact with each other. Created 
at this interaction are a shock stronger than either of the original two, 
a contact surface, and a weak expansion wave. As soon as the sonic Mach 
number is reached at the nozzle throat, no more expansion waves can pass 
through, and therefore only an expansion fan of finite width can travel 
to the left through the supply tube. A short time after the flow in the 
throat chokes, a Q shock is formed just downstream of the nozzle throat 
and slowly travels through the remainder of the nozzle and the test sec- 
tion. A steady flow is established behind this Q shock, and the start 
time is therefore defined as the time when the Q shock passes out of the 
test section. Some of the Q and P characteristic curves are traced in the 
wave diagram presented in Figure 2. It is seen that the paths of the P 
characteristic curves are fairly uniform throughout the wave diagram, 
but the paths of the Q characteristics do not follow any systematic 
patterns in the early stages of the flow. After establishment of steady 
flow, both characteristics follow a uniform pattern; i.e., at each x sta- 
tion the slopes of both characteristics are constant in time after estab- 
lishment of steady flow. At any time, we can recognize the flow regime 
at any position in the facility by observing the directions of the Q 
characteristics. Subsonic Mach numbers are characterized by the Q 
characteristic curves traveling to the left, sonic when traveling in 
the vertical direction, and supersonic when traveling to the right. 

Figures 3 through 6 present the wave diagrams for test Mach numbers 
1.4, 1.7, 3.3, and 3.0. The wave diagrams for all supersonic cases are 
similar to Figure 2, and therefore the Q and P characteristic net is not 
shown in Figures 3 through 6, because to enable one to draw the charac- 
teristic curves requires a print-out at each time step by the computer 
program. This print-out takes much more computer time and results in 
volumes of paper. 

Figure 7 plots the static pressure as a function of time at a sta- 
tion in the test section for the supersonic range of test section Mach 
numbers calculated. The drop in pressure to a constant value after the 
shock passes this station thereby establishes the starting time defined 
as the beginning of the period of steady flow. 

For the subsonic cases, the flow is choked downstream of the test 
section by choking flaps which are approximated by assumed boundaries 
that account for the flaps and blockage due to the model support and 
diaphragm cutter mechanism. Cases were calculated in which subsonic 
and sonic Mach numbers of .315, .7, and 1.0 resulted in the test section 
during the period of steady flow conditions. For the subsonic, as in 
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the supersonic cases, only an expansion fan of finite width can pass 
through the throat since no more expansion waves can pass through the 
throat after choking occurs. Steady flow is established in the test 
section only after the highly nonlinear effects of area changes on the 
slopes of the P and Q characteristic variables at each x station 
upstream of the throat have died out with time. This effect is illus- 
trated in Figure 8, which shows the pressure as a function of time at 
an axial station in the test section. The starting times for the sub- 
sonic and sonic cases are defined as the times when steady state con- 
ditions are established throughout the test section. 

Starting times for all of the cases that were calculated are plotted 
in figure 9 as a function of the test Mach number. It is seen that the 
subsonic and sonic Mach numbers require a longer start time, since the 
steadying of the flow properties in the test section is a gradual pro- 
cess whereby, for the supersonic cases, a strong shock passes through 
the test section that instantly establishes a steady flow behind it. 

The starting times for the supersonic cases depend, however, on the 
speed at which the shock travels through the nozzle and test section. 

This speed is governed by the local flow conditions upstream of the 
shock and the strength of the shock. 

Figure 10 presents the calculated Reynolds number per foot in the 
test section during the period of steady flow for the range of test 
Mach numbers and for the maximum charge pressure of 48.64 atm assumed 
for the facility. It is seen that the highest test Reynolds number of 
approximately 2 x 10 s per foot for the facility will occur in the vicinity 
of a test Mach number of 1.4. Figure 10 shows one case which was cal- 
culated for a test Mach number of 2.04 in which the charge pressure was 
assumed to be 18.03 atm. This reduction in charge pressure resulted in 
approximately a linear reduction in test Reynolds number for this 
particular Mach number. The effect on starting time by this reduction 
in charge pressure was found to be negligible, however. 

The static pressure in the test section during the period of steady 
flow is presented as a function of Mach number in Figure 11. The expected 
trend of decreasing pressure with increasing Mach number is shown. 

To check the effects of the settling chamber on the start time, one 
case for M^g^ = 2.04 was calculated in which the settling chamber was 
removed. A comparison of these results with the corresponding case with 
a settling chamber showed a negligible difference in the start time. 

The preliminary analytical results that have been calculated by 
this procedure will be compared with experimental results when the facil- 
ity is completed and measured data made available. Future analytical 
investigations to include the effects of boundary layer and mass removal 
of gas through the walls in the test section for transonic Mach numbers 
are planned. 
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AERODYNAMIC BOUNDARIES OF FACILITY WITH THE MACH 2 NOZZLE 
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AERODYNAMIC BOUNDARIES OF FACILITY WITH THE MACH 1 NOZZLE 
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FIG. 7. PRESSURE AS A FUNCTION OF TIME 
AT A STATION IN THE TEST SECTION 
FOR SUPERSONIC TEST MACH NUMBERS 
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START TIME AS A FUNCTION OF TEST MACH NUMBER 
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APPENDIX A 

Derivation of the Characteristic Equations 


The basic differential equations (1) through (4) can be written in 
nondimens ional form (see Definition of Symbols) as: 

Continuity 


1. 5 (PA) , 1_ S(PuA) , 
A St A Sx 


* = 


0 . 


(A-l) 


Momentum 


Su , Su 
St U Sx 


af S(ln P) 
7 Sx 


(A-2) 


Equation of State 


a 


2 



T. 


First Law of Thermodynamics 


S 


S 


o 




— In p # 

7 


(A-3) 


(A- 4) 


Some useful relations that result from the combining of equations (A-3) 
and (A-4) are 


p/p ? = e 7l7-ms-s 0 ) 


(A- 5a) 


p = a 2 7/ (7-1) e -7(S-S 0 ) 


(A- 5b) 


2/(y-l) -7(S-S 0 ) 

p = a ' e 


(A-5c) 


55 



Due to the relations (A-5a) through (A-5c), the continuity equation 
(A-l) can be expressed as 


— 2-r + -2-r u & + a & - - au - a 20"*! 

7 - 1 at 7 - 1 ox dx dx St 


+ ^d + “S)-^ <r3)/(r ' 1> e7(s ' So) - 


and the momentum equation as 


du du , 2 Sa 

St U Sx 7 - 1 a Sx 


a2 S + £ - 


Adding equation (A-7) to equation (A-6) results in the following 




. +a (7-3)/(7-l) e 7(S-S 0 ) + f> 


where 


P = 


- 1 


a + u. 


The substantial derivative in the x,t-plane is defined as 


D _ S , _S_ 

Dt St U Sx 


Let the differential operator 


+ _ S , s S 

— = ^ + (u + a) ^ 


(A-6) 

(A-7) 

equation: 

(A- 8) 

(A- 9) 

(A- 10) 
(A- 11) 
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represent the change of a parameter in the direction 


dx 

3T u + a 


in the x,t-plane. Equation (A- 8 ) can then be written as 


6 + P _ ... 3(ln A) . S(ln A) . , , . DS , 8 + S 

sE au ' a -^dT- + a <7 - X) D? + a bT 


- \|ra 


(7-3)/ (/-I) e 7 (S-S 0 ) + f> 


Subtracting equation (A-7) from (A- 6 ) yields 


$ + (u . a) | . . au ^ . a Mail + ra (k . i 


at 


Sx 


at 


Dt 7 ax 


. ^(r-3)/(r-D e 7(s-s 0 ) . £> 


where 


Q = 


7 " 1 

Let the differential operator 


a - u. 


c) , s 5 

m aF + <u ' a) s 


represent the change of a parameter in a direction 


dx 

dt 


= u - a 


in the x,t-plane* Equation (A-14) can then be written as 


(A- 12) 

(A- 13) 

(A- 14) 
(A- 15) 

(A- 16) 
(A- 17) 
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5_Q 

5t~ 


= M ln _ A). . a d(in 4) , 

aU Sx St + 3(7 


n DS + 5 - S 

1} ^ + a iT 


^(7-3)/ (7-1) e 7(S-S 0 ) 


f. 


(A- 18) 


A third relation is required to complete the system for the three 
dependent variables a, u, and S. Given by the entropy condition which 
must be prescribed for any problem, this relation is given here in a 
general form as 


nq 

— = F(a,u,S,x,t). 


(A- 19) 


Equations (A-13), (A-18), and (A-19) form the system of equations that 
must be solved for the three dependent variables u, a, and S. These 
equations are solved by a step-by-step procedure along the characteristic 
curves P, Q, and S in the x,t-plane. The details of solving these equa- 
tions are presented in other sections of this report. 
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APPENDIX B 


The Fortran Program and Its Input Data 


1. Preparation of the Input Data for the Computer Program 

The computer program and the subroutines employed by the pro- 
gram are written in Fortran IV language for the CDC 3200 computer in 
Part 2 of this Appendix. To run the program, the following input cards 
must be prepared: 

Card No. Columns Fortran Symbol Format Description 

1 1-14 TOL E14.8 Tolerance for which all 

solutions calculated by 
an iteration procedure is 
satis fied . 


15-28 TOLQ 


29-42 P2P1 


43-56 A2 


E14.8 Maximum difference 

between the values of the 
characteristic Q for neigh- 
boring Q waves (see Sec- 
tion X, Paragraph B) . 

E14.8 Ratio of the pressure on the 

left-hand side of the dia- 
phragm to that onthe right 
before rupture (p^/p*). 

E14.8 Nondimens ional speed of 

sound on left side of dia- 
phragm before rupture; see 
equation (33). 


2 1-14 G 


15-28 XREF 


E14.8 Ratio of specific heat ( 7 ), 

The program assumes that the 
gases on both sides of the 
diaphragm have the same 
value for 7 . 

E14.8 Total facility length (L 0 ) 

by which all distances are 
nond imens ional ized . 
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Card No. 

Columns 

Fortran Symbol 

Format 

Description 


29-42 

XD 

E14.8 

Diaphragm location (x^) . 


43-56 

XTEST 

E14.8 

Axial location (x^g^) in 
test section for which 
pressure and Mach number 
are determined and printed 
out. 


57-59 

M 

13 

Number of cards used in 
describing tunnel geometry 
(see Section X, Paragraph A) 

3 to M+2 

1-13 

XE 

E14.8 

Axial station (x^ for which 
diameter and curve parameter 
are given. 


15-28 

DE 

E14.8 

Diameter (d^) of the axially 
symmetric duct at x* 


29-42 

EPS 

E14.8 

Jo 

Curve fit parameter (e“) at 
xg (see Section X, Paragraph 
A). 

M+3 

1-14 

DT 

E14.8 

Nondimens ional time interval 


15-17 

NUMT 

13 

Number time steps to be 
located . 


18-20 

NPRNT 

13 

Printout frequency. 

M+4 

1-15 

T 

E15.8 

Nondimens ional time . 


For T > 0, the program incorporates a restart procedure that 
requires information from the last calculated time in order to resume 
calculations. This information is contained in a number of cards that 
have been automatically prepared at the termination of the previous run. 
Therefore, for T > 0, the input data, beginning with card #M+4 plus the 
remaining required input, is automatically prepared at the termination 
of the previous run. 

For T = 0, card #M+4 is the last input data card. 
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2. The Fortran Program auu uuuiuuimea 


PROGRAM DEC* 3 

C CONST ANT TIME INTERVAL METHOD 

dimension z r< to) ,cki ) .xq( i ) .oq(i) 

COMMON G,UT, fOL . tolo 

COMMON XE(25),nt<25).EPS<25),YC(?5.3>,M,XTL,XTR,Xl) 

COMMON XXG< 350 ) .PPG < 350) ,QQG< 350 >, AAG( 350) ,uUG< 3511 > ,SSG( 350 > 
COMMON XG(350 ) , pg<350 ) » QG ( 350 ) , AG(350 ) ,UG(350 ) . SG ( 350 ) 

COMMON Z1 (350 ) ,/2( 350 ) ,XG8(350 ) , XG&(350 ) 

COMMON PM1.PW1.PPM1.FPW1 

COMMON PVPI, Al, J1,P1.C1,S1,A2,U2,P2,Q2.S2, A3 , U3 , P3 , Q3 , S3 
COMMON NNG< HI ) ,OOW(10) .QOM(10 ) . NO ( It) ) .UW(10 > ,QM(10 ) ,NNQ r.NOT 
COMMON MMCS ( 1 0 > , MMCST, MCS( 10 ) , MCST 
COMMON XSHOCK(IO) 

COMMON X CONTACT! 10),* CON TACT <10 ) 

PEA U (60, 88) rCL.TOLO.F2Pl.A2 

88 FORMAT (4E14.8) 

RE A 0(6 C, 90)G . XREF , XU. XTEST »M 
90 FORMAT ( A E 1 4 . 8 , 13) 
v TI-3T = X T EST/x rfF 
C = 2 . * G / ( G - 1 , ) 

_PRNT=1 

C G IS THE RATIO OF' SPECIFIC HEATS 

XPEF" IS THE REFERENCE LENGTH WITH DIMENSIONS (TOTAL LENGTH) 

XD IS THE DIAPHRAGM LOCATION WITH DIMENSIONS 
M IS TH= TOTAL NUMBER OF GEOMETRIC POINTS READ IN 
WR I TE ( 61 ,69 ) 

89 FORMAT ( //3X , 29HC0NSTANT TIME INTERVAL METHOD/) 

WR I T E ( 6 1 ,95) X REF 

95 FORMAT(3x,2HLr,El5.8.lX.16H(REF6RENcE LENGTH) ; 
y D = x D-/ x R E F 
NR I TE ( 61 , 94 ) x IJ 

9 4. FORMAT < /3X,5hxn/ 1. = .ET 5.8, IX, 19H( DIAPHRAGM SiATI ON )////) 

W R I T E ( 6 l , 3 0 D ) 

300 FORMAT ( 34X.15HT0N.NEL GEOMETRY/) 

WRI T E ( 6 1. , 3 0 1 ) 

3d 1 FORMAT ( 4X.1HI ,6X. 5HXE ( I ) , 7X, 10HXF( 1 ) / XREF , 8 X , 5HDE i I ; , 10X, 6HEFS( I ) 
1 6 x , 1 1 HE PS ( I )/XPfcF/) 
no HO I - 1 , M 

P E A 0 { 6 0 , 1 on ) XE r I ) , DE ( I ) , EPS ( I ) 

100 FORMAT ( 3E14 . 8 ) 
vK D =xE( I ) 

EPSP=EPS< I ) 
ve ( i ) = x c < I ) /x.xr-F 
FPS ( I ) =EES( I ) /XREF 

an WRT r.-( 61 . 302 M , XEP. XE ( I ) . DE( I ) , 6PSP, EpS( I ) 

302 FUHMAI ( 3x , I3.5F15.8) 

C VE ( I ) IS VALLE OF X (WITH DIMENSIONS) AT STATION I 

0 TE ( 1 ) IS DIAMETER (WITH DIMENSIONS) AT Xfc ( I ? 

FPS(I) IS SMALL INCREMENT ON EACH side of Xfc ( I ) FOR WHICH THE 
C CUBIC CURVE FIT FOR THE QIAMETER IS USED 

v T i = X E ( 2 ) - E P S ( 2 ) 

DO 3l5 1=2, M 

■ IFmtST-XK i ) ) 3 1 6 . 3 1 6 . 3 1 5 

315 CONTINUE 

316 rrEsrs=i 

I T L S f L = 1 - 1 

TR1=ITESTR+i 

L'O 398 IR=!R1,M 


61 



IFCUEC I R)-I)fc( IR-l) >356.398,399 

398 CONTINUE 

399 yTR=(XE(If<)+XE( l P ■ 1 ) > / 2 . 

F-X IMAX=CX 6 < l I eS I R)-XE< 1T6STL) > /20 . 

TxOMAXs . 015 

nxIMIN=0XlMAX/3. 

rxOMlN=DXOMAX/2. 

DOI MAXsTOLO- 
PO I M I N= TOLO/3 . 

- rOOMA**.l 6 

rOUrtlNsDUOMAX/?. 

tall geomcm, xe, oe, bps, rc> 

u.RlTfc(61,303)- 

303 FORM A T-L / / / 5 X > 53HC0EFF I C I ENT S FOR PARAriOL 1 C CURVE FIT BETWEEN SECTI 
l.r ns > 

vRf TEC61.304) 

30 4 FORMAT ( 3X , 19PD = C1 + C2‘S I *C3*S 1 **2 ) 

WRI (6(61,30?) 

.50 5 FORMAT (3X,66ESI=<X- (>El I ) - EPS ( .1 ) > >/XREF AND XE(I>-EE<1> .LE. X .L.E 
1 . Xt < I ) +FE ( P/> 

WRIT EC 61 , 306 ) 

306 FORMA T ox , 1H 1 . 5X . 2HC1 , 13 X , 2HC2 , 13 X. 2H<;3 ) 

R 1 •- M - 1 

no 81 1 =,, mi 

81 UR ! TEC 61. 307)1. YC(I,1),YC(I, 2) , YCC 1 , 3) 

>-•,07 FORMA T(3x. 13.4615,8) 

Al = l . 

1 ■ 1 =0 - 

17 = 0 . 

p! = ?.*Al/(0-l. 

61=0. 

P2=2,*A2/CG-l . ) +U2 
n2=2.*A2/(G-l. )-u2 

REA-BC 60 , V 8 ViLT , NUMT , NFRNT 

98 FORMAT (cl4.8,2I3) 

C FT IS THE TIME INTERVAL 

c mjmt is the number of t waves to compute 

C NPRN'f IS the MUMPER FF T WAVES BETWEEN PRIM OUT 

u r rr e < 6 i , 2 oo ) r 

?oo format < //-Sx. fhg= , 615 .8. 2X , 7h< gamma ) ) 
ukl rb( 6 i,?oi) tolo 

201 FORMAT (3X,5hT0L0=.E15 . 8 ) 

UR 1 T E ( 61 , 202 ) P2 3 l 

202 FORMAT! 3x.6hP2/P1=.El5. 8 ) 

UR 1 TE ( 61 ,203 )DT 

2 0 3 . FORMAT ( ix ,.3*4C T = , El 5 . g ) 

WRI TE(61 ,11 0 ) 

110 FORMAT ( 1H1.) 

READ ( 60 , 86 ) T 

86 FORMAT ( El 5 , 8 ) 

I E ( T - T t J i_ >1, 1 . 3 
1 ruNflNut 
^ =5u 

<r u-sNUMBRR E(> I NTS 06SIWE0 IN INITIAL EXPANSION FAN 

! PGO=0 

U: :'j { - 0 
ECS I =1 
M1E-N 

C ALL T WAVF t ( N Q E ) 

■ J .4® 2 - - 

T = T + 1) T 

FAIL T WAVE 2 (NQE.NCS .N,JJ) 

PCS ( 1 ) =NC$ 
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. K--( 0|P*Ai4^1><U , 4 2 , 4 3 

4? jRRNTsI) 

4. 3 FONT t- N4 i 

0 0 TO 3? 

J awj-mis 

PEAU(6 n, ? 0 ?>)lpi?0. NOT, MOST, N 

?fl 3 F-04<MA| < 1 1 . ? I 3 > - -- 

F E A IJ < 6 0 , 2 0 6 ) (MCS( I ) . 1=1. MOST) 

ES AU 0 64- ( M t) ( T->-,-T-=4-rNai4 — 

206 FORMAT < 10 I 3 > 

P0--7- T-= l , NOT 

7 FEALi(6'),2Q7)Cw( I ) 

? q 7 r^iputAT(pj p;_a) . . . 

TO 4 KQ = 1 . N 

p 6 A 0 < 6 ,Wj-0-M4-rX-6-(-K^T V -R6 (-K-&-W-4G LK-64 -SG !-*&T 

106 FORMAT ( 13.4615.3) 

MT-pK-iTT-.3*-(-PG(KG)-0G(KaW ----- 

AG(Ki)) = (R-l. )*<RG(<Q)+OG<KQ) )/4. 

4 -CONTINUE 

F2 = S0 ( 1 > 

IK !PG0)P,8.Si 

-G.EAXU 6 G-. 1.0.7 4F M 1 *.R ty 1 _ 

107 F 0 * M A T ( 2 F 1 5 . 2 ) 

9 JNUE - _ ... - ...... ... ... 

„'PRi\iT = 0 

3 2 FONT tNllE . 

TO 29 wr=l,NLMT 

jRRN f = JPRNT+1 

jfXG(D-=XG C 1 ) +OT *AUG U ) -.A£( 1) ) ... . 

TFlXXGU) >500.500, 510 

*v ii (i ru. 6ui i.=juju 

301 7 T ( l ) = 2 . * X G ( I )/0T + UG( I ) - A G ( I ) 

_AAG ( 1> = AG ( 1 ) ._ 

302 AGUESS = A AG ( 1 ) 

CiX.U.=AAGXJLl- . ... 

CALL lNTER(2.in,l,ZT .XG.C1.X0) 

CALL -J-NTER-C2^-N , 1 , XG . CG . XQ. QQ ) - «.... 

AAG(l) = (fi-l. )*OOU )/2. 

.11- (4.8.S (AAG(l)-AGOESS)-T OL )5Qi. 503, i?02 

303 vy.G(l)=0. 

PPU(1)=2.*AAG<3 )/(G-l. ) . . 

COG ( 1 ) = PPG ( 1 ) 

I 1 I I f: f 1 ^ = fi _ . . 

SSGlUsS? 
v«tt( 1 )= 0 . 

yGQM ) = x to ( J ) 

-P4±— ^ o4 J = 1 , 1 0 

504 7T( I ) = X3( I )+CT*(UG( I )-4G< I ) )/2. 

-++£{ -2->-*AG44>~ -- - . - - 

MJG(2)=JG(1) 

3 0 3 AGUEOS = A AG ( -24 — 

1 GUESS = JUG ( 2 i 

— -X-X4< 2>-XG( 1 ) ♦ i>T * ( UG ( 1 ) +AGI1 4+UUGT 2 ) + AA&4.2 J )/2 . 

Cl(l)=XXG(2)-DT*(UOG(2)-AAG(?))/2. 

CALL IMER<2.10.1, 7T,XG,C1,XQ> - - 

CALL I N T P R ( 2 • 1 0 , 1 » X G » C) G , X Q » QQ ) 

UUG-U2-) -=4-4641 >~04-UtA>-/-2-, 

AAG<2)-(G-1.I*(PG<1)+G3(1))/4. 

4 F +AUS < JyG<2 i-uGtlE-SS ) - TGU5G6,50-6»505- 

50 6 IF( AMS( AAG( 2 )-AGUESS)-TOL >507,50 7,505 

-^7-.-PA4A(-2-)=2,-*-AAG-(-2 )/T-G--lf4*UUG-(-2 ) 

CQC,< 2) =2 . *AAG (2 ) / (3-1 . >-UUG(2) 

F -S P4 2 4 *SG4 -14 



V(5h(2>sX6(l ) 

-vs-a-(y-)-=x(n i ) 

ro 508 1=1, N 

-4F4 Xii< -I HIOO > >-5nB,508,509 

508 COOT I N U t 



.10 = 2 

. -a*)- T O- 5 11 - 

510 CUNTiMUr 



C 0 G t 1 )=Q2 

I IJG ( 1 ) = J? 

. SS6U->=S2- 

VG*<1 )=0. . 

- - - *G*±< 1 > = XG( 1 ) 

IQ - 1 

iCQ^V- - - 

511 r ONTINUt 

.. no 2 UhN 

7 1 ( i ) = X G ( I > -»• L T * ( L) G ( I ) + AG( 1 ) )/2. 

3 - 7 - 2 < i >rx-GM KGH >*0T/?.- 
II « = 1 

o0~l 
„JO = l 
jC^ = l 
-2JC;s=i 

400 CONTINUE 

10= 1 0^4 - — — ■ . 

M=I» 

K^KQ + 1 

CAI.L PnlNT0( HO. KKQ, T , IQ,KQ, N. 0 , DTP ) 

PI XsABSf XXG< lO)-XXG(Uj-l)) 
r I O = A B S ( U 0 0 ( iQ)-OQG(IQ-l)) 

IF <XX(i( 1 fa ) - X X G f I 0-1 ) >427,427,401 

401 IF ( XXG( 10)-XTL)4D2,406,406 

A\X? j-F (-01 X-#X0M4X > 44U3 ,443-, 4 1 7 

403 iFMJIG-OnOMAX >404,404,417 

404 1 F ( 0IX-QX0M1 N) 40 5, 4 05, 4 22 

4 0 5 1 f ( D I Q - Ut-OM I K ) 4 1 8 , 4 1 8 , 422 

406 TF< x <&< U?) .S 1 . • XTL . ANC . XXG( IQ ) . LE . XTP MO 7, 411 

407 lF([HX-UylMAX)4 0 8, 406,417 

* 0 8 - T-P T-U- PO~0 0 OP A X M-0 9 , 4 0-9-, 417 

409 TF(DIX-DXIM1K)410, 41 0,422 
4-10 IF ( DIG- DQ I MI M 41 8, 416 ,4 22 
411. TF(XXG(U‘)-1.)4)2,412,412 

412 IF( DI X-DXOMA 0416, 416, 41 3 

413 KKX=IQ-1 

_F A — FT) i KT A 0-C H-O-v N > - - 

ro 414 XKJrKM , I 0 

IF < x X G ( <K J ) - 1 . >414,415,415 

414 CONTINUE 

- 415 lUrKKJ 

A M = I Ci 

4-1 6- — f-Ak -b— -S-U-P F ND ( I N-,-1 P GO ) 

GO TO 6] 

417 r.Al.I. P-0 IMADC < I Q • N > 

M\t - 1 0 

- .... - GO TO 422 

418 JF(KQ-NQ(J0) Ml 9, 4? 3, 4 [9 

40J9 lEAALh-MZSUiC *4 14 2GU 42 5-^420 

*♦2 0 IF MQ-NM21 ,23,421 

4.2JU .1-0=10-1 

GO TO 4 u (j 

— 4-2 2 I F(KU-WOUfj) M24.4 23.424 . 
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423 Tall OOCK(T,IQ,i<Q,J0,JJQ. jcs,jjcs,n,nn. IPGO' 

JxZA Li (KQ-ri-C4-XJC£ > .>.4.24 > 49 5— A26-. - — - 

"25 CALL CONT(T, fQ.KO, JQ. JJQ, JCS,JJCS,N. IIO.KKQ.NN, IPGO) 

426 _,.IF(AU-N)4U0.23.400 ■ 

427 CALL CR0SS< IC.KG. JO.wJQ.NEWQ) 

IF(NteWU)422, 422,428 — ■ ■ 

426 jPPNT = N°RNT 

23 TF< IPG0)26.24 ,26 

2.4 f ALL PHUCK( IC,N,.MN > - 

IF(XXG(NN)-1. )5 1.51,63 

83 I F ( LLii ( I (j ) » AAG ( I 0 ) lB4- t -£A+25 ....._ . . - 

84 N N = N N - 1 

F -Ai . l . S u 8fc.QL l< M 4- r A) - ) - - • 

TpGO = l 

Rfl T O. 31 

25 IQ=NN 

. CALI SaJife-NDt-i-CU (CQ, iUMi,. I P 00. l - 

26 CONTINUE 

J F t 1 PGO ) 51 . 51 ,431 

4 V TF(XXG(NN>-1. >50,51.51 

5 0 CALL Suit NO ( K N » M ) 

51 CONTINUE 

- MslilT =JJ3-1 

PMCST = JJ(,S-1 

iFLOfiRNT-OPOM >28,27,23 

27 FAI L PRINTOUT (T.NN. IFGO.XTEST) 

- ^ ? 9. N T = U 

25 CON I I NUc 

FALL MOv/PtNN, N, IP-GO) - 

CO TO (4?9,29),SSUTCPF(4> 

2 9 FONflNU s 

"29 uRlTr<61.30fl>T 

3nS FOEMAT(3x,5HTlMfc=»irl5,3//) - 

WRl TE(61 ,101 > 

-- l‘)l FOPMAT(4PX, IPX ,4-5 X , 1FP ,14x,lWQ, 1 4X.lHA-.4-4x, 1-HU, 14X , IPS, 13X.4P MACH, 
1 1 0 X , 6HPPPSS ) 

AHACH = I!3( 1 ) / A a ( I ) 

3 0 WRITE: (61, 105) I,XG( I > ,*PG( I ) , QG< I) , AG< I ) ,UG( I > ,SG( I ) ,AP ACH , PRESS 

1Q5- FORMAT < 3X , 13,«E16.8) - - - 

UR I TP (62,06 )T 

-VR-f-T-E ( 6 2-r 2 0 9 ) -FPG O ) -ROT-,-MOOO,-0 - - 

URi TE(62,206) < MCS( I > » I =1 » MCST ) 

uO-i4P4-62-T 2-0-6) -4«K44 t J-4T4G T) - 

HO 34 I =1 , NQ r 

34 w«I TE<62,207 )QW< I ) - 

TO 35 1=1, N 

3-9 — u - R I T P ( 6 2 , 1 0 -o )-pi -X-G ( 1 ) , P G ( -I-)-,-QG(4--)-,-S G ( I ) — 

I F ( I PGO ) 37 , 36 ,37 

36 L.R l et ( 62 . 10 7 ) PMl .P wl- - — - 

■37 CONTINUE 

PNG 
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' cuWROUYlNlEi PC iNTADDt IQ.N) 

r i mens i ON xBUl.A Bt D.URM > , CB ( .1 1_JL C i 1 ) . UC 1).,.SC(1) , X S ( l ) . SS ( 1 1 . Cl . . 
irl),C2(l).C3 l l>.Z3<300)»AC(l) 

COmilM-JL-D T . I OLiYULQ — . - 

rOMMON XE(25).D6(25).6PS(25), YC(25»3),H,XT .XTR.XD 

common xx G( ?5q i_._p.pg (.ISO J_, J10GX33-aj S£G ( 35DJL .. .. - 

rOMMON Xfi(350 ),PG( 350 ),QG( 350 ).AG(350).Ul*( 60),SG<350) 

rOMMON 71(350i./2( 350 1,.X.GB(55 1U_^G.QL55.fl) _ 

ro m 1=1 . n 

.1-0. j3 ( Ll = XG(1 1 (J.fiC I ) b AG1L) ltDT ... _ — — - 

CXX=XXG(!Q)-XxG( IQ-1) 

RAOXs-l A AG ( I £)- A AG4 1 Q-=l > U-QXX — 

rufiX = (iJUG< IG)-UUG( I 0 - 1 > ) / DX X 

i F ( x X G ( 10) -GE.XT L-A NE.XXG t.i m . I..E » X Tr 1 1 1 2 

1 TUI. OC=.15> 

2 T0I.QC = T01Q 

3 rn = Q QG( lCll.C.Q&l IQ.-1Y. .... — 

ro 11 1=2,4 

r[)Q=iio/Ai 

JF LABS ( Ujjflj ,6t • TQLQ.C111 ,12 

11 CONTINUE: 

CO TO 13 

12_I±L=_L=J 

1.3 i-=nxx/Ai 

V X G ( 1 A ) =XXG < iQ) 

PP lii LA 1-P1GXXQ-) — — 

CQG ( I A ) =QQG ( 1 Q ) 

AAG (I A ) = AAG( LXP. 

I. UGC I A ) =l)UG ( * Q ) 

— . <?SLU LAJ-SSSG-(-l-O) — - - 

yGH< 1A)=VGR( • Q ) 

XGU.C 1A ! =XGQOQ) - 

I Q = 10 — 3 

10= 1 0+1 

XA = XXG( ICl-D +H - - ----- 

CALL AREA(XA’DAA) 

AA=AAG(10-l)+OADX*(XA=XXG(lq-l)) 

I: A - 0 1 j G ( IQ -1 )*niJL>X *±M s x 1G LL£Lt JJ 1 

14 aGUESS=AA 

jC= JC+1 

GJJJ..1 = VA^, 5* ^UA+ A_A ) *1 T - 

C3(l)=XA-.5*llJA-AA)*CT 

r? ( 1 )=XA-.s*Oa* DT 

CALL INTER <2>N.1.Z1»XG. Cl, XB> 

T F ( XB ( 1 ) -X G ( 1 ) >20.20.21 — 

20 A B ( 1 ) = A G ( 1 ) 

i b a ) =u g a > — 

C B ( 1 ) = S G ( 1 ) 

CO TO 22 

21 CONTINUE 

r Ai I 1 NTER ( 2 • N , 1 . XGj_AJL._XJL,.ABJ - 

CALL I N TE-R ( 2 ’N.l.XG.LG.XB.UB) 

C Ai l lMTFRf7.-N.l -XG . S G . XB . SB-1 
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?? CONTINUE 

- CALL INTER (2*N.1,Z-2.XG.C2,XS> 
l£XXSI.LL=JL GX i. ) „) .23 . 2 3. .2.4 


23....SS(1.IlS-SJSa 1) — 

CO TO 25 

—~3A CONTINUE — - - — 

CALL INTER<2.N.1.XG,$G,XS,SS) 

2- 5 - COAU -Uats 

CALL INTER(2»N.1.Z3,XG.C3,XC) 

- CALL I N T-fcR ( 2-> N -1 ».XG ,-A G .. X-C , ACT . 

CALL INTER(2'N.1.XG.IG.XC,UC) 

CA-L-L -4N T -fe it( 

CA=SS( 1 ) 


-C-AJJ — AREAJxatj ) ,DAHA- 


CAI.L AREA (XC ' ,DAC) 

P P = 2- . * A acl-M- Lfi^-L^U Ufe U.4- 

rc=P.*AC(i)/iG-i.)-uc<D 


-SB ( 1 ) U-2~ 


CA=QC+(-UA*AA.UaA-UC( 1)*AC(1)*DAC)*DT/2.+( A*AC(l))*(SA-SCa))/2. 
AAS-tl».:;.l.-l.»fP**QA)/4 


I A=.5*(PA-GA ) 


15 lRITE(61,102> 

hRITF(61,103 MlGUESS.L A, AGUESS, AA 

10T rnRrtAT<3X.7HUfillPSR=.F15.B.?X.3HUAs.Flg.H.g .7HAGUPRSs.PL5-fli2X.3HA 


1 A= , tl5 . D/ ) 

-JLQ_m_.l.SL_— 


16 CONTINUE 


17 TF( ABS( AA-AC-ESSI-TOL >18,18, 14 

18 Y XG ( I Q ) = X A 

PPG ( IQ)=?.*AV(G-1. )*UA 

COG ( IQ)=2.*AA/fG- 1 . )-UA 

AA6( IQ)=AA 


-U UG..C -LQ L=TJ A. 
PSG( IO)=SA 


yGB ( IU) sXb(l l 

X GO ( IQ)=XC(1' 

1 9-f-CALU-NOS 

IO=I A 

ftfcl UR N — 

END 
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3200 FORTRAN ( 2 . 1 . 0 ) / ( RTS > 


SUBROUTINE RO»NTQ< I IG.KKQ.T, IQ.KO.N, ID.DTP) 

DIMENSION C1<A).XB<2),AB<1),UB<1)»SB<1>»C2<1),XS<2>,SS<1>,DTT<2) 
DIMENSION USS ( 1 ) , PB8 ( 1 ) 

COMMON a.DT.TOL.TOLO 

COMMON XE(?5)»DE(25),fcPS(25),YC(25»3)»H,XTL.XTR,XD 

COMMON XXG(35U) ,PPG(350 ) .OQG(350 ) . AAG(350) »UUG(350 ) ,SSG(350 ) 

COMMON X G ( 3 5 0 J ,PG(350 ) . Q G ( 3 5 0 ) , A G ( 3 5 0 ) , U G ( 3 5 0 ) » S G ( 3 5 0 ) 

COMMON Z1 (350 ) . Z2<350 > . XG8<350 > . XGO<350 ) 

COMMON RMl.PWl.PPMt.FPWl 

COMMON P2P1.A1.J1, PI, 01, SI, A2.U2»P2.Q2»S2. A3,U3,P3,U3,S3 
COMMON NNOflO > , QOWdO ) ,QOM(10 ) ,NO(10 ) ,QW(10 ) ,0M(1D ) , NNQT, NOT 
COMMON MMCS tl'J),MMCST,MCS(10),MCST 

uc=n 

CAl L AREA (XG ( *0 ) , DAQ ) 

A A = A G < K Q ) 

IIA = UG ( KO ) 

19 AGI'ESS = AA 
IIGU6SS = UA 
JC=JC+1 

XA=XG(K3)+nT*lUA-AA+lG(KQ)-AG(KQ))/2. 

C1(1)=XA-.5*(UA+AA)*CT 

CALL iNTER(2,N,i,Zl.XG,Cl.XR) 

I F < ID)?3,20,2P 

20 IF ( XR( 1 )-XG(in21 ,21, 22 

21 A 9 ( 1 ) = A 2 
UB ( 1 > =M2 
SB(i)=S2 
GO 10 28 

2 2 nTP=OT 

CALL INTER(2,N.1,XG,AG,XB,A8) 

CAl L INTER(2.N,1, XG,FG,X9,PBB) 

UB(1)=P98(1)-^.*AB(1)/(G-1. ) 

CALL INTER(2,N.1,XG,SG,XB,SB) 

GO TU 2 8 

25 OXX = XXG( I IQ)-XO(KKO) 

IF ( XB(1 )-XG( Kl\Q) ) 24 , 2 4 , 22 

24 OTT (l) = OAaxX**( I IQ) ) / (UA + AA-DXX/DT ) 

IF( XA-XXGC 1 10 ) )1, 1, 25 

1 S A = S G ( K 0 ) 

XR(1)=0. 

GO TO 41 

25 I • B ( 1 )=MJG( I IO)-OTT(1)*(UUG( I IQ)-UG(KKQ) )/DT 
AB(1)=AAG(II0>-0TT(1)*(AAG<IIQ)-AG(KKQ))/DT 

PTT(2) = (XA.xxu(IIQ))/((UA4-AA + uB(1)*AB(1))/2.-DXX/DT) 

IF ( DTT ( 2 ) )59 , 05 , 58 
55 IF (OTT (2)-DT)05,65,59 
59 XB( 1 ) =XA-L)T* ( UA + AA ) 

6n CALL INTER(2.N,1,XG,AG,XB,AB) 

CALL INTFR(2» | 5.1,XG,LG,XB,U8) 

XB ( 2 ) =X A -DT * (UaVaa + UE (l)+AB(l))/2. 

IF ( XB( ?) -XG(KNO) ) 61,61,62 
6 1 XB( 1 )=X3(KKC ) 

A8 ( 1 ) = A G ( KKQ ) 

UB(1)=U3(KKC) 

SB ( 1 ) =SG ( KKG ) 

DTP=DI 
GO TO 28 
69 CONTINUE 

IF ( A8S( X B ( 2 ) - * 8 ( 1 ) )-T0L)64,64,63 
65 XBC 1 )=X3(2) 
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GO TO 60 


64 m P = D T 

CALL INTFR<2 ,N,1. XG,SG. XB,SH> 

GO TO ?S 

65 CONTINUE 

IF( A3S(DTT(2 )-PTT(l ) )s .1F»07)27,?7*26 
26 DTT(1)=DT7(2> 

GO TO ?9 

2 7 S B ( 1 ) = S S G ( I IQ>-DTT(1>*(SSG( I I Q ) -SG ( KKQ > ) /DT 

XB( 1)=XXG< no>-DTT(l)*DXX/DT 
DTP=DTT(1 ) 

25 PBs2.*A3(l)/(ti-l. >+UE<l> 

C2<1 >=XA- .5*D I * U A 

CALL 1 N-T ER ( l , >M . 1 , l? , X 3 , C2 , XS ) 

I F ( 10)99,29,34 
29 TF(XS(i ) - XG ( 1 > >30,30*31 
3n S S f 1 ) s S 2 
GO TO 37 

3 1 CALL INTEH<2,«M,XG,SG,XS,$S> 

GO TO 37 

39 IF ( XS< 1 ) -XG(K»'Q ) >33,33. 29 
33 rm(l)=<xA-xx*(ilG>)/(UA-DXX/DT> 

3 4 LlS = *JuGf n 0 > - D * T ( 1 ) * ( L OG ( I IG)-UG(KKQ) ) /DT 

DTT(2> = (XA-xX^(II0))/((UA-*.uS>/2.-DXX/DT> 

IF(UTT (2) ) 5 1 » P 1 . 5 0 
5 n IF(UTT(2)-DT>37,57,51 
5 1 X S f 1 ) - X A - U A * D I 
59 CALL INTEM2,‘'M,XG,AG,XS,USS) 

XS(2)=XA-DT*(OA+USS(l > )/2, 

IF (XS(?)-XG(KftQ) ) 53 » 5 3 » 5 4 
53 SS(1)=SG(KKG) 

GO TO 37 

5 4 IF ( ABS < XP< 2>-XS <.1 ) >-TOL> 56,56.55 

55 XS(1) =Y S(2) 

GO TO 52 

56 CALL I N I E R ( 2 , N , 1 , X G , S G , X S , S S ) 

GO TO 37 

57 CONTINUE 

IF <A8S(DTT( 2) -DTT(t) )=. 1E= 07) 36,36*35 
35 DTT(1)=DTT(2) 

GO TO 34 

3 6 SStl )=SRG( ! !0>-DTT<l)*(SSG( I I Q > -SG< KKQ > ) /DT 
37 SA=S3 ( 1 ) 

CALL AREA(XA,UAA) 

CALL ARtA(X0( J-) »DAB) 

PA=Prt+( -UA*AA»DAA-LJB(i )*AB( 1 ) * D A B ) * DTP/2 . + < A A + ABU ))*(SA-SB( 1 > )/2. 
OA=QG(KQ ) + ( -lici( KO) +AG ( KQ > * D AO- U A* A A *D A A ) *DT / 2 . + < AA + AG (KQ) ) *< SA-SG< 
1 K 0 > ) / 2 . 

A As (G-l , ) * ( P A + Q A ) / A , 

1 1 A = . 5 * f PA -QA ) 

TF ( J C - 2 9 ) 3 9 , 3 9 , 38 
3 5 WRITE (61,100) 

3.00 F 0RMATr/3X, 33<*1T0LERAI^CE CANNOT BE MET IN POINTO) 

LID JFF =UGUESS-OA 
AD1FF=A3LESS-*A 
WRlTfc<61,'l01)AA,UDIFF, AD IFF 

11)1 FORMAT { 3X,2HX = ,E15.8,2X,6HUDIFF = # E15.8*2X,6HADIFF=,615.8/ ) 

GO TO 41 

3 9 IF ( AdS( JA-UGU«=SS)-T0L > 40 * 40.19 
4 n IF ( AHSf AA-AGIJCSS)-T0L >41,41.19 


41 X X G ( I Q ) = X A 
I ! l.i li ( I 0> =MA 
A AG ( I Q ) rM 

PPG( 10) =2 . *A A / ( G-l . ) +UA 
OOG( IQ) =? . *AA/ (G-l . ) - U A 
SSG ( 10) =SA 
XGR( IQ)=Xb(l ) 

XG0( 10 ) =xG( «Q7 

RETURN 

END 
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320 0 FORTRAN < 2 . 1 . 0 ) / ( RTS ) 


SUBROUTINE: T *A VE 1(NQE> 

COMMON 3,UT, TUL . TOLO 

COMMON Xfc ( ?5 ) * HE (25 ) , EPS ( 25 ) , YC ( 25 , 3 ) , M, XTl . XTR, XD 

COMMON XX6(350) .PPG(350 ) #QOG(350) . AAG(350 ) »UUG<350 ) ,SSG(350 ) 

COMMON XG(350> ,PG( 350 ) , QG< 350 ) , AG( 350 ) »UG(350 ) ,SG( 350 ) 

COMMON 21 ( 35 0 ) , 7.2(350 ) , XGB< 35 0 ) , XGQI350 ) 

COMMON PMl.PVU, PPM1.FPW1 

COMMON P2H1,A1,U1»P1»Q1»S1»A2*U2»P2»Q2»S2*A3»U3»P3»Q3»S3 
P1P1=1. 

C = 2.*G/(G-1. ) 

S2 = Sl-AL0G(P2Rl/A2**C ) /G 
P1=2.*A1/(G-1 . )+Ul 
01 = 2. *A1 /(G-l* ) - U 1 
P?=2.*A2/(G-1* ) +U2 
Q2 = 2 . *A2/ (G-l • )-U2 
S 4 = S 2 
P4 = P2 

C ITERATION PROCEDURE EEGINS 

1)3 = 1 .2 
1 I.I3G-U3 

As - ( G+t . )*(ll-U3)/Al 
7=( A + SORT ( A**2 + i6 . ) ) /4 . 

A3A1=S0RT((2. + (G-1.)*Z**2)*(2.*G*Z**2-(G«1.)))/(<G+1.)*Z) 

AG1=2.*B*Z**2/(G+1.)-(G-1.)/(G+1.) 

AG?=( ( 1 . + ( G-l • ) * 7* * 2/ 2 . )/( ( G + l. ) * Z * *2/2 . ) )**G 
DS = ALO(5(AGl*AP2)/(G*(G-l. )) 

P3P1=A3A1**<2.*G/(G-1,))*EXP(-G*DS) 

P4P1=P3P1 

A4=(P4P1*EXP(u*(S 4-S1 )))**( (G-l . ) /(2. *G) )*A1 
i,I4 = P4-?.*A4/(u-1. ) 

U3=U4 1 

IF(A8S(U3-U3G>-T0L)2,2,1 
P CONTINUE 

C ITERATION COMPLETED 

P M 1 = 7 

04 = 2 . *A4/ (0-1 . ) -U4 

P4=2.*A4/(G~1* ) +U4 

A3=A3A1*A1 

03 = 2. *43/ (G-l- >-U3 

P3=2 . *A3/ (G-l . ) +U3 

s3=si+ns 

PW1=U1+A1*PM1 
T = 0 . 

WRITE(61,115 ) I 

115 F0RMAT(//3X,5HTIM£s,E15.8//) 

WR I TE ( 61 ,117) 

117 FO»MAT( 5X , 1M I » 7X, 1HX, 14X , 1HP, 14X, 1HQ, 14X.1HA , 14X,1HU,14X, 1HS/ ) 
A Q = N Q E - 1 
DO= ( 02-04 ) / A Q 
DO 3 KO=l,NCE 
01 sKO-1 
XG ( K 0 ) = X 1) 

PG ( KQ ) = = 2 

OG(KQ)sQ2-DO*'4I 

UG(KQ) = .5*<PG(KQ)-QG(K0) > 

AG(K0) = (G-1. )*(PG(KO)+OG(KO) )/4 . 

SG(KQ)=S? 

WRI TE(61«101 )rtQ» XG(KC ) > P G ( K Q ) , Q G ( K 0 ) , A G ( K Q ) * U G ( K Q ) »SG(KQ) 

101 F0RMAT(3X, I3.0E15.8) 

3 CONTINUE 
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WR I T E ( 6 1 .131 ) 


131 F0RMAT(/3X.42"PR0PERTIES ON LEFT SIDE OF CONTACT SURFACE) 

WRITE (61, 132 ) R4.04.UA, A 4, S4 

13? FORMAT (3X .?HP=,6l5.8,2HQ=,El5.8,?HU=,E15.8.2HA=,El5.e,2HS=,Elt>. 8/) 
WRITE(61,133) 

133 FORMAT* /3X.43HPR0PERT IES ON RIGHT SIDE OF CONTACT SURFACE) 

WRITE (61 .132)P3,03.U3,A3,S3 

WR 1 TE ( 61 , 13 A ) 

134 FORMAT ( /3X.37HPR0PERT IES ON LEFT SIDE OF SHOCK WAVE) 

WRI re ( 61. 132) ^3, 03. U3, A3. S3 
WRITt(61.135)^Ml,PWl 

135 F()RMAT(/3X,2nnp SHOCK MACH \IUMBER= . E15 . 8 , /3 X , 17HP SHOCK VELOCITY*. 
1F15 , 8 ) 

WRITE (61. 136) 

136 F0RMAT(/3X,38MPR0PERTIES ON RIGHT SIDE OF SHOCK WAVE) 
WRITE(61,132)H1,01.U1.A1.S1 

WR! TE ( 61 ,110) 

110 FORMA |(1H1) 

RETURN 

END 
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3200 FORTRAN ( 2 . 1 . 0 ) / ( RTS ) 


SUBROUTINE T WAVE 2 ( N UE , NCS . N . J J ) 

COMMON G.DT.TUL.TOLO 

COMMON XE<25) » D6 ( 25 > , EPS < 25 > , YCI25.3) , M , XTL . XTR , XD 

COMMON XXG(35U) , PPG ( 3 50 ) , QQG ( 350 ) » A AG ( 350 ) » UUG < 35u ) » SSG ( 350 ) 

COMMON XG(350)»PG(350)»QG(350)#AG(350)*UG(350.)»SG(350) 

COMMON 21 ( 3*50 A , 72(350). XGB(350),XGO(350) 

COMMON »Ml,PWl,PPMl.FPWl 

COMMON P2P1 ,A1. U1,R1,Q1,S1. A2,U2,P?,02*S2» A3»U3,P3,Q3,S3 
TC= JJ-1 
T = TC»OT 

XXG(1)=XD + DT*H;2-A2> 

PPG( 1 ) =p? 

OQG< 1 ) =02 
A A G ( 1 ) = A ? 

OUR ( 1 ) = J? 

SSMl )=S? 

DO 3d (=?,NuE 

V X G ( I ) = X D + D T * M J G ( 1 ) - A G ( I ) ) 

PPG ( 1 )=PG< I ) 

COG ( I ) = QG< I ) 

A A R ( I ) =AG( I ) 

IHlR ( I > = UG( I ) 

3D SSG ( I ) = SG ( I ) 

XCS=XU+U3*DT 
NO? = 9 
AQ? = NQ2 

P X ? = ( XCS-XXG (NOE ) ) / AC 2 

K 1 = N O E ♦ 1 

NCS=NQF+NQ? 

no 31 K=K1,NCP 

XXG(N)=XXG(K-A)+DX2 

PPG ( < ) =a>pG(nqc) 

OQG ( K ) = 0GG ( \ Oc ) 

A AG ( K ) = AAG(NQc ) 

UliR ( K ) = UbG ( N O c 1 
SSG < K ) = SSG ( NOc 1 
3 1 CONTINUE 

xSP = XU + UT*Pvil 
N O 3 = 1 0 
A Q3 = j NO 3-1 

PV':( XSR-XXG ( NCS ) ) / AG 3 
K2=NCS+1 
N = N C S + N G 3 
DO 3? k=k2,\ 

AK-K-K? 

xxh(K)=XXG(NC-J)+DX3*AK 
PPG ( K ) = p 3 
0(JG ( K ) =03 
A A n ( a ) = A 3 

iiug < K ) = J 3 

3 ~> SSG ( K ) = S3 

PPM1=PM1 
PPW1=PW1 
OR 1 1 E ( 61 . 115 ) i 

11 o FORMAT ( //3X , 5"T I ME=,E15. 8, IX, 3H (.)//) 

WR I TE< 61 , 117 ) 

1.17 F0RMAT(5X.1HI»7x.1HX.1*»X,1HP,14X,1HQ.14x,1HA,i4X.1HU,14X,1HS/) 
PO 33 1=1. NCS 

33 WR1TEC61, 101)1, XXG(I)»PPG(I) »QUG( I ) , A AG ( ! ) ,UUG( I ) .SSG ( I ) 

101 FORMAT (3X, 13. °E15.S) 
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WR I TE < 61 # 137 ) 

137 FORMAT(3X.15HwO\ITACT SURFACE) 

MCSR = i\lCS*l 
no 34 I = N C S H » N 

3a WR 1 Tr < 61. 101 ) * , XXG< I ) .PPG ( I ) , OOGf I ) » A AG ( I ) » UUG< 1 > . SSG ( I ) 
WR?TE(61.13l3)PPMl,PPl»l 

13=> FORMA r< /3X.201P SHOCK HACH NUMBER*. E15 . 8 , /3X, 17 HP SHCCK VELOCITY 
1E15. 8) 

WR r TE ( 61 . 110 ) 

11 n FORMA fflHl) 

no 36 1=1, N 

XG( I ) = XXT-( I ) 

PG ( 1 ) =PP(.i ( 1 ) 

OG< I )=OOG< I ) 

A G ( I >=AAG( i ) 

UG ( 1 >=MJG< I ) 

3 =5 SG < I ) = S $ G ( [ ) 

PM1 =PHM1 
PW1 =PPW1 
RETURN 
END 



3200 FORTRAN ( 2 . 1 , 0 ) / ( RTS ) 

SURHOUTINt CO*T(T, IQ.KO, JQ. JJQ, JCS, JJCS.N, I IO,KKQ,NN, IPGO) 
common a.UT.TUL.TOLO 

COMMON XF (?5 ) »f)£(25 ) , bPS(25) , YC(?5,3 ) , M, XTL, XTR,XD 

COMMON XXb(3 I 5 u ' , PPG ( 2 £ 0 ) ,QGG(350 ) , A AG < 350 > » UUG ( 3bQ ) , SSG(350 ) 

COMMON Xtt(350i ,PG( 35 0 ) ,00(350 ) , AG<- 350 ) ,UG( 350 ) ,SG( 350 ) 

COMMON 71 <35 0 J , 22( 350 ) ,.XGB( 35 0 ) , XGQT350 ) 

COMMON 3 Ml ,Pwl, PPM1, FPW1 

COMMON P?Pl.Al.ul,Pl,Gl,Sl,A2,U2,P2,Q2,S2,A3,U3,P3,Q3,.S3 
COMMON MivO( 10 > , QQW(10 ) , OQK( 10 ) , NO( 10 ) , QW(10 ) , QM(10 ) , N NOT, NOT 
COMMON MMCS( 1»J ) . MMCST , MCS( 10 ) , MCST 
COMMON XSHOCKMO) 

COMMON XCUNT ACT (10 ) , (.CONTACT ( 10 ) 

1 CALL CONTACT ( *G,KQ, JC, JJQ, JCS,JJCS,N. I IQ,KKQ,NN, IPGO, T) 

K = 10—1 

IF (KO-M) 16, 13* 1 3 
10 IF ( KO-MCS ( JCS ) ) 2, 1 , 2 
P I(J= 10 + 1 
K 0 ■= K 0 + 1 
N N = J 0 

CALI. PHI NTO ( I i Q , KKQ, T , I Q , KO , N , 1 , OTP ) 

IF(KQ-MCS( JCST ) 19,1,19 

19 I F ( X X G ( i 0 ) - 1 . M , 3 , 3 

X CALL SU?H:ND< lu,KQ,N,NN, IPGO) 

GO TO 13 

a I F ( K 0 - m 3 ( j Q ) ) o , 5 , 8 
5 I F ( X X G ( 1 Q ) • X X « ( IQ-1) >20,20,21 

20 T Q = I 0 - 1 

21 CALL QSmOCK( I«,KO, JO, JJO, JCS, JJCS,N, I 10, KKQ) 

IF(XSHOCK( J*.0-1 )-XCOKTACT( JJCS.-1) )6,6,7 

A CALL CSCHOSSOtT, IQ.KC, JQ, JJQ, JCS, JJCS, I 10, KKQ, N) 

7 00 TO ? 

B I F ( K 0 - N ) 9 , 1 3 , i 3 
9 TF(OTP)ll ,11, AU 

10 IF ( XXb < \ Q ) .XX<J ( I 0-1 ) ) 11 , 11 , 12 

11 T Q = l Q - 1 
00 10 ? 

IP IF ( AHS(UTP.UT)-T0L)13, 13,2 
IX CON I I N'Jfc 
R H r u R N 
F N f • 
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3200 FORTRAN ( 2 . 1 . 0 ) / ( RTS ) 


SUBROUTINE: CONTACT ( I C , KQ , JO , J JO , JCS , J JCS . N , I IQ, KKQ, NN, IPGO, T) 
DIMENSION XP(<d ,AP(1),UP(1),SP(1),C1(1) ,XB(2) ,UB(1) , A B ( 1 ) , SB(1) , XC 
l(2).UC(l).AC(A).C2d),XSd)»SSd),SCd) 

DIMENSION x6NU(1),AEND(1),UEND<1),SEND<1> 

COMMON 3,DT,TUl,TOLQ 

COMMON XF(25).DE(25),EPS<25) , YC(P5»3).M.,XTL,XTR,XD 

COMMON XXG(35« ) . PPG ( 350 ) , QQGC 350 ) , A AG (350 ) , UUG< 35 0 ) , SSG( 350 ) 

COMMON XG(350 dPG(350),QG(350 )»AG(350)»UG(350),SG(3'50) 

COMMON Z3 (350>,Z2(350> , XG6 ( 350 ) , XGG ( 350 > 

COMMON RMl.PWi.PPMl.PPWl 

COMMON P2P1, Al.Ul.Pl.Ql.Sl, A2.U2,P2,Q2,S2, A3 , U3 , P3 , Q3 , S3 
COMMON NNQ(10J,QQW(10),QOM(10),NO(10),OW{10),QM(1U),NNQT,NOT 
COMMON MMCS(1U ) , MMCST, MCS(10 > , MCST 
COMMON XSHOCKdO) 

COMMON XCONTAOT(10),LCONTACT(10) 

NCS=MCS ( JCS ) 

NCS2=MCS( JCS+1) 

CALL QDI.S(JQ,XDISL,XCISR) 

SL=SG(NCS) 

SR = SG ( NCS+3 ) 

01=fcXP( ( G- 1 . >* (SL-SRJ/2. ) 

F=1 . + fcX 3 ( (G-l ♦ )*(SR*SL>/2, ) 

DTO=OT 
IIL = UG(NCS) 

AL=AG(NCR) 

JC = 0 
1 IJ = UL 
A = A L 
JC = JC + 1. 

IZ = 0 
IP = 0 
UR = UL 
AR=AL/01 

XR=XG(NCS)+CT*(UR+UG(NCS+l))/2. 

XL = X R 

XP(i)=XR-OT*(UR-AR) 

IF ( XR- 1 . 12,32*32 
? continue 

I F ( XP ( 1 ) -1 . ) 3 A , 32 , 32 
31 CONTINUE 

CALL INTFR(2,N, 1, XG, AG, XP, AP) 

CALL InTER(2,N,1,XG,L(J,XP,UP) 

XP(2)=XR-DT*(UR-AR + UFd)-APd))/2. 

IF( JCS-MCST)33.36,36 
3=5 IF (XP(?)-XG(NLS2) ) 36 , 82, 82 
36 CONTINUE 

IF (ABS(XP(2)-XPd ) )-TUL)4,4,3 
3 VP (1 ) =X P ( 2 ) 

GO TO 2 
32 CONTINUE 
VEND ( 1 ) =1 . 

CALL 1NTER(2,1Q,3.,XXG,AAG,X6ND,AEND) 

CALL I.NTER(2»JQ»1,XXG,UUG,XEND»UEND) 

CALL INTER(2, 1 Q,l, XXG.SSG.XEND, SEND) 

NN= I U+l 
KQ = N 
I PGO= 1 

XXG ( NN ) =1 , 0 0 0 U 01 
AAG(NN)=AENC(1) 

UUG(NN)=UEND(1) 
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SS 6 <NN)=SfcNCU) 


PPG(NN)=?.*AAP(NN)/<C- = 1. )+UUG(NN) 

Our, ( NN ) = 2 . * A Au ( NN ) / ( G si . ) ■U'.JG ( NN ) 

XGHf NN)=() . 

XGQ(NN)=n . 

GO TO 7 5 

4 IF< jCS-4CST)8.i.,83,83 
8 1 NCS2=NCS( JCS + D 

IF(XP(l)-XG(NwiS2) >83,82.82 
8? UP(1)=UG(NCS2> 

AP(1)=A(J<NCS2> 

SP(1)=GG(NCS2> 

SMT = lJG ( M0S2 ) 

SM?=(UR-AR*LP(l)-AP(l ) )/2, 

TO=T+( v^-XG(NPS2)-DT*SMl)/(SMl-SM2) 

XU = XW + S'1?*( TU-T ) 

DTOsT-TQ 
XP( 1 >=X0 
1P = 1 

GO t 0 6 

8 3 C All. Ii»TFR(2,W.l,xG,Si;,XP,SP) 

IF(XP(1)-XDIS*>6,5.5 
3 K1. =NU(J3) 

1.7 = 1 

XP(1)=X3(K1) 

AP ( 1 )=AG(K1) 
i IP ( 1 ) =I)G(K1 ) 

SP( 1 ) =S3 i K1 ) 

A 0P=2.*AP(l>/(u-l. )-UP(l> 

CA! L AREA ( XB , PAR ) 

CAM.. APEa(XPU) .DAP) 

OR=QP+(-HP(l)«AP(l)*tAP-UR*AR*DAR)*DTQ/2.+< AP ( 1 ) +AR ) * ( SR-SP ( 1 ) )/2. 
Cl(l)=X!.-DT*(UL + AL>/2, 

GAIL lNTFRt2,N.l,Zl.XG,Cl,X8) 

JF(XB(l)-XDlSl-)60,60,64 
6(1 Jl = NNQ< JJQ.D + 1 
JI=.N0( J3-D+1 
OUH=( 01.13 ( JL >-UG( JI ) )/DT 
nAR=(AA3(JL)-«G(JI))/l>T 
DXO=(QOH( JJK-D+OW( JG=1 ) )/2. 

X8(l >=XL+(UL +a L )*(XXG( jl)=xl)/(ul + al-dxq> 

61 TBT=(XP(l)..XX'a( JL) )/CXO 
UB(1)=UJG< JL)”DUB*TBT 
AB(1)=AAG(JL) + 0AB*TBT 
nxPW=(Ui.*AL + UO(l)+AB(l) )/2. 

XB ( 2 ) =XL + DXPW*< XXG( JL )-XL)/(DXPW = DxQ) 

IF ( A B S ( XB(2)-XB(1) )-T()L)63, 63,62 
6 P XB(1)=XB(2) 

GO TO 61 

63 SB ( 1 ) =S5G ( JL ) + ( SSG ( JL ) -SG ( J 1 ) ) * T B T / D T 
OTP=-THT 

GO TO 65 

64 CONTINUE 

IF ( JJCS-2J89 , ° 4 , 84 

84 NCG1=MCS( JCS-A)+1 

I F ( XB ( 1 )-XG(NOSl> >85,05,89 

85 MCCS1=H1CS( jJPS-D+l 

PXX = XC0MTACT(PJCS-1>-XG(.NCS1 ) 
nTTl=(XL-XCGNiACT(JJCS-l))/(UL+AL-DXX/DT> 

86 UB ( 1 > — 1 1 C 0 N T AC • ( JJCS-1 ) - DTT1 * ( UCONT ACT ( J JCS - 1 ) -UG ( NCS1 > )/DT 
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A8(l >=AAG(NCC^1)-DTT1*( AAG<NCCS1)-AG(NCS1) )/DT 
DTT2=(XL-XCCNl ACT( JJCS-l) ) / < <UL + Al«*UB<l)*AB< 1 > >/2, -DXX/UT) 

\F ( AHSf DTT2-DI T1>-.1 Er07>88,88, 87 
87 1)TT1 = 0TT2 
GO TO 86 

8r sen. )=s3( ncsi > 

Xb ( 1 ) = XCONTAC ( (JJCS-l ) -DTT1*DXX/DT 
OTPri)TTl 
GO TO 65 
86 CONTINUE 

CALL INTFR(3,*,1,XG,IG,XB,UB) 

CALL. INIE-R(2,'^1,XG, A(i,XB, AQ) 

CALL INTER<2, N,1,XG,SG,XB,S8) 
r>TP = OT 

6 S P8 = 2.*A3Cl)/(ta-i. )>UB(1) 

CALL AREA(XB(A),DAB) 

PL=PB+(-Ub(l)*A8(l)*CAB-UL*AL*DAR)*DTP/2,+(AL*AB(l) > * ( SL-SB ( 1 )) /2 . 
AL=(G~1 . >*<PL + 0R)/(2.*E) 

UL=PL-2.*AL / ( «-l. > 

IF (JC-30)2l, 2^.20 
2H UDIFFsu-UL 
AD!FF=A-AL 
WRI TE<61 , 90 ) 

9n F0PMAT(/3X,28HT0LERAM;E NOT MET IN CONTACT) 

WRITE* *1,91) XL, AD1FF,UD1FF 

9i F0PMAT(3y,3hXL=,El5.e,2X,6HADIFFe,Bl5.8,2X,6HUDIFFs,ei5,8) 
go ron 
2 1 continue 

TF ( AbS(UL-L')- I 01)9,9,10 
6 IF(ABS(A! -A)-‘0D11, 11,10 
In GO TO 1 

li continue 
K'CSL 2 I 0 + 1 
NCSR=NC5L+1 
HMCS( JJCS)=NCPL 
VCONTACT* JJCS > = XR 
U CONTACT (JJCSI=UR 
XXG(NCSl)=XL 
PPG(NCSl)=PL 
AAG(NCSl)=AL 
UUG(NCSlI=UL 

0QG<NCSl)=2.*Al/(G- 1 . ) -UL 
GSGfNCSL )=SL 
XG«(NCSL.)=X6(D 
XGO( NCGu )=0 . 

XXG(NCSR)=XS 
QQG<NCSR)=OR 
A AG(NCSR)=AR 
UUG(MCSR)=UR 

PPG(NCSR)=2. *AR/(G-1. )+UR 
SSG(NCSR)=SK 
XGP(MCSR)=0 . 

XGO(NCSR)=XP(A) 

I QrNCSR 
NN= IQ 
I 10= IQ 

TFf 17)51, 51, 5U 
5n KQ=Kt~l 
K K 0 = K Q 
GO TO 52 

5l IF( IP)54,54,5»> 

53 KQrNCS? 


GO TO 52 

54 no 6 6 I=1,M 

IF ( X G ( I ) -XP( 1 n 66, 66 . 67 

66 CONTINl.'c 
KQ = N 

K K 0 = K Q 
GO TO 52 

67 K 0= I - 1 
KKOsNCS+1 

5P JCS=JCS+? 

JJCS=JJCS+1 
75 CONTINUE 
return 

ENO 


3200 FORTRAN (2 . 1 . 0 ) / (RTS ) 


SUBROUTINE COG«(T, IQ, KQ, JQ , J JQ, JCS , J JCS , N,NN, IPGO) 

COMMON' 3, DT, TU|_, TOLQ 

COMMON Xfc<25>#DE(25),fePS<?5),YC<25.3>,M,XTL,XTR,XD 

COMMON XXG< 35J ) , PPG< 3 50 > , GQG( 350 ) , AAG(350 ) ,UUG(35(j > ,SSG(35d ) 

COMMON KG ( 35 0 > . PG( 350 ) , QG( 350 ) , AG ( 350 ) , UG( 350 ) , SG( 350 ) 

COMMON i ! ( 35 0 > , 12 < 350 ) , XGH ( 350 ) , XGQ ( 350 ) 

COMMON »Ml , PWi , PPM1 , FPW1 

COMMON P2Pl,Ai,<Jl,Pl,Ql,Sl,A2,U2,P2»Q2»S?>A3>U3,P3»Q3,S3 
COMMON Mi*(l< 10 ) , QOW( 10 ) , QOM( 10 ) . NQ 1 10 > , 0W < 10 > , QM < 10 > , NNQT-, NOT 
COMMON MMCS(iy)»MMCST,MCSC10),MCST 
COMMON XSHOCK'10) 

COMMON X CONTACT ( IQ ) , L CONTACT ( 10 ) 

17 CALL GSHOCK< N.KG, JQ, JJQ, JCS, JJCS,N. I IQ,KKQ) 

N M = I Q 

1F( JJQ-2)20,2U,18 

1R IF ( XSHOGK ( JvQ-3 ) -XSHCCK( JJO- 2 ) ) 19, 19 . 2 0 

19 CALL UCRUSSC ( ' , I 0, KO, JQ, JJO, .JCS, JJCS . I I Q, KK(J, N) 

20 CONTINUE 

21 10=10+1 
KQ = K(J + 1 
NN= 1 0 

CA! I.. PC INTO ( I i Q , KKQ, T , I 0, KQ, N, 1 , DTP ) 

T K ( X X G ( I U ) - 1 . >23,22.22 
2? C A • L Sii a t-ND ( I W , KQ, N, NN , I PGO > 

GO TO 25 

23 I F ( KfJ-NU ( JO ) ) 45 , 24 , 25 

24 GO TO 17 


25 

If 7 

(KO-mCS< JCS> ) 2 6 , 

,10,26 

10 

IF 

( X X G ( I ( i ) - X X « ( 10- 

■1) >11,11,28 

1.1 

IQ 

= 10-1 



GO 

TO 2 3 


25 

IF 

(KO-N) 27,28*28 



27 CONTINUE 

IF(0TP)1B, 15, 14 


14 I F ( X X G ( I Q ) « X X U ( IQ-1) ) 1 5 , 1 5 , 1 6 

15 I 0 = ! U - 1 
GO TO 21. 

If, CONTINUE 

TF(AriS(i)TP«CT>-TOL >26,28,21 

2 R CO N TlNHe- 
PETilRiM 
F N 0 
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3?nO FORTRAN ( 2 . 1 . 0 ) / < RTS ) 


SUBROUTINE CS rt OCK ( IQ.KQ, JO, JJQ, JCS, JJCS, N, I IQ.KKO ) 

DIMENSION XL< M » AL<l>*UL(i> ,SL<1> ,XC<2>, AC<1> #UC(1> ,SC<1> 

COMMON fi,DT,TUL,TOLQ 

COMMON XF(?5)'D£(25) , EPS (25) ,YC<25.3> , M , XTL . XTR , XD 

COMMON XXG<35<J ) .PPG (350 ) ,QQG(350 ) , AAG(350 ) ,UUG(35U ) , SSGT 350 ) 

COMMON XG(350>,PG(350)»QG(350)»AG(350)»UG(350 ) ,SG 1.350 ) 

COMMON 31(350) , Z2(350)»XGB(350),XGO(350) 

COMMON - > Ml,Fwl,PPMl»FPWl 

COMMON P2P1 .A1,U1,P1,G1,S1, A2.U2,P2'Q2'S2» A3,U3,P3,03,S3 
COMMON NNQ(10),OOW(10),GOM(10)»NO(10)»QW<10),GM<10)»NNQT.NOT 
COMMON MMCS(IO) .MMCST,MCS(10),MCST 
COMMON XSHOCK(IO) 

COMMON XCONTAUT(10 ) . L CONT A C T ( 1 0 ) 

NU=NO( JQ) 

call xdiscomjq, jcs.xuisr.nqocs) 
jc=o 

GQSWsQV ( JQ) 

10 GUf SS=OQSW 
17 = 0 
JC=JC+1 

HA = ( OQSR+QWC Jw> )/2. 

V L ( 1 )=XG(NL )+‘)T*WA 
XR = XL ( 1 ) 

I F ( JJCS-i ) 2 » 2 * 1 

1 IF (XCONTACTt JJCS-D-xLd) )2,4,4 
? IF( JJQ- 1)5,5, P 

3 IF(XSHOCK( Jw.Q-l)-XL(l ) >5,4,4 

4 AL<1)=AG(NL) 
l'L(l)=l.'G(NL ) 

SL ( 1. ) =SG ( NL > 

GO TO 6 

=5 CONTINUE 

CAI L INTER(2, 10.1,XXG,AAG # XL,AL) 

CALL INTER(2, 1Q,1,XXG«UUG,XL'UL) 

CALL InTER( 2, *0,1, XXC- , SSG , XL, SL ) 
f CONTINUE 

PL = 2.*Ai_(l)/(^-l. )+UL(l) 

QL = 2.*Al(1)/(«*-1. )-UL <1 ) 

PQSM = -( aOSR-ULcl) )/AL (1) 

0 U = 2 .*( 1 . -OGSM* *2 ) / ( ( G + l. )*QQSM) 

UR = UL ( 1 ) +D(j* At- ( 1 ) 

A R = ( S Q R T ( < 2 . + ( G*1 , )*C«SM**2)*(2.*G*Q0SM**2-(G-1, )))/((G + l.)*CQSM)) 
1 * A l_ ( 1 ) 

XC( 1 )=XR-DT* (UR-AR) 

11 CONTINUE 

CALL INTER(2,N,1,XG, AG.XC, AC) 

CALL INTER(2,N.i,XG,LG,XC,UC) 

XC(2)=XR-DT*(UR-AR*UC(1)-AC(1) )/2. 

I F ( A3S(XC(?)-*C(1) )-TOL >13,13.12 
12 X C ( 1 ) = X C ( 2 ) 

GO TO 11 

13 CALL INTER(2,N,i,xG,SG,XC.SC) 

IF (XC(l)-XDI SR >17,17,14 
Id 1F(NQ0CS-1)17»15,16 
13 Kl=NO(J3+l) 

17 = 1 

XC(1)=XG(K1) 

AC(1)=AG(K1) 

UC ( 1 ) =U3 ( K1 ) 

GC(1 )=SG(KJ ) 
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c 0 H» i 7 


If Ki=MCS<J(S> 

17 = 1 

XC<l)=Vi5(M ) 

a t: ( i ) = a s < K 1 ) 

UCML ) =HK«K1 ) 

SC ( l ) =sS(ki ) 

17 CALL AP-A (XC ( 1) . DAC) 

CALL AREA(YR,uaR) 

0C=?. ^AC(l)/(«-l . )-UC(l) 

OR =2 . *AR/(G-1 * ) -UR 

net =(UR-nu/A*-a> 

CALL SHOl.M CE*-, 'JOSM, LR AT ( ARAT.PRAT.DS) 

SR-bLM )+DS 

OR=OC+r - aC( 1 ) *UC<1 )*C AC-AR«tJR*DAR )*DT/2.*(AR+AC(l))*(SR-SC(l)>/2. 
DEI = ( QR-PL ) / A L ( 1 ) 

CAt L bHUUK(CEt- f QOS^#LRAT,ARAT,PRAT,DS) 

OOSrtsUI. C 1 ).AL ll ) *QQS^ 

!F< JC-40 )19. IP, 18 
IP PIFF=AHS(OOS^-GJESS) 

WRTfE(6i»9Q) 

9n FORMA U /4X, ^"TOLERANCE CANNOT BE MET IN GSHOCK) 

WR l Tb( *1 , 91 HUFF 

9i FORMA r f3X,1 AH * OLERANCfc MET=.Fl5-8/> 

CO T n ?n 
IP CONTINUE 

IK(AtfSOOSM-GUESS)-TCL) 20 .20.10 

20 OQw( JjO)=UQ$H 

DOM JJO ) =OOSM 
no 23 1=1 , I C 

J F f xxG f I >-XR >*1 ,22 . 22 

21 CONTINUE 
2? 10=1 

NMO ( JJO ) = IQ 
X SHOCK ( JJQ) =XL ( 1 ) 

y x < I Q ) = x L ( l ) 

Ppf- < 10 ) =PL 
000 ( 10 ) = QL 
A A r. ( I Q ) = A L ( 1 ) 

I'UC.C 10 ) = i 1 L ( l ) 

ssc( lQ>=vin ) 
x G h ( I Q ) = l . 
y go ( I Q ) = n . 

10=10*1 
I io=iu 

Y X r, ( 1 0 ) = X R 
OQM IQ) =CK 

A A G < I0)=AAG(T«J-1)*ARAT 

HUG ( I Q ) =2 . *AA^ I Q ) / ( G = 1 . ) »OQG ( 10 ) 

PPU{ lO)=?.*AAb( lQ)/(Gil. )*UUGf IQ) 

SSG (JQ) =SSC( I 0 - 1 ) +DS 
XGM ( IQ) =(- . 

YfJO ( 10) = XU(1 > 

I F ( f 7 ) ? 8 , 2 A , 2 ^ 

2 7 K Q = K 1 - 1 
K K Q = K 0 
go to 2 p 
2 B CO [ ' T I N u E 

Du 2 3 1=1 . N 

T F ( X G ( I ) - X C < 1 >>23,23.24 
23 CONTINUE 


KiJ = N -1 
KKO = KCj 
GO TO 

24 K Q - I - 3 

K K 0 = M Q ( J 0 > + 1 

25 J J O = J j O + 1 

ja= jq*i 

RETURN 
F NO 
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3200 FORTRAN ( 2 . 1 . 0 ) / ( R T S ) 


SUBROUTINE SUFEND( IQ.KQ,N,NN, IPGO) 

DIMENSION XENU< 1 ) , PEKU( 1 ) , OEND(l ) »SEND(1> 

COMMON G.DT.TUL.TOLQ 

COMMON XE ( 25 ) » DE ( 25 ) , EPS ( 25) , YC( 25 , 3 > » M , XTL . XTR, XD 

COMMON X X G ( 3 5 U ) » PPG (3 50 ) , QQG(350 } » A AG ( 350 ) » UUGC35 0 ) , SSG(350 ) 

COMMON XG(350>»PG(350)»QG(350)»AG(35Q)#UG(350),SG(350) 

COMMON 21(350) ,Z2(350 ) , XGB ( 350 ) , XGQ ( 350 ) 

X E N O ( 1 ) =1 . - 

CALL I NTER ( 2 » * G . 1 , XXG , PPG , XEND , PEND ) 

CALL INTER<2. *0.1, XXG, QQG, XEND. QEND) 

CALL I NTER ( 2 . *0.1, XXG, SSG, XEND, SEND) 
nn= i g 

X X G ( N N ) = X E M C ( 1 ) 

PPG ( NN ) =PENC ( ■*■ ) 

OQG ( NN ) = Q E N C ( a ) 

AAG(NN) = (G-1. )*(PPG(NN)-*-QOG(NN) )/4. 

UUG(NN) = .5*(PPG(NN)-CQG(NN) ) 

SSG(NN)=SEMC(i) 

XGP ( NN ) = 0 . 

XGQ ( NN ) =0 . 

KQ = N 
I PGO = l 
RETURN 
END 
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3200 FORTRAN < 2 . 1 . 0 ) / ( RTS ) 


SU H R(HJT ) Nfc SUUEND(NN.N) 

COMMON S.DT.TOL.TOLQ 

COMMON XP ( 2s ) • DE ( 25 ) .EPS (25), YC(25>3) .M,XTL,XTR,XU 

COMMON XXli(35» ) .PPG (350 ) ,QQG(350 ) . AA6(350 ) . UOG( 360 ) ,SSG( 350 ) 

COMMON XG(35 0 ) , PR (350 > .QG<350 ) . AG (350 ) ,UG(350 ) . SG<350 ) 

COMMON Z1 (350) , *2(350 ), XGB< 350 > .XGO(350 ) 

DIMENSION X8U).AB(1),UB(1).SB(1).SIB(2) 

JC« = 0 
N N = N N + 1 

AM2= ( XXG(NN-1 > -XG(N) > /DT 
XXI' ( NN ) =1 . 

CALL APEA ( XXGI NN > , DAE ) 

HUG C NN ) =UUG ( N‘H-l ) 
m UfiUfcSSsUUGf GNi 
20 AM1 = IJUG(NN) 

SIB(1) = <3. .■AM1*UT-XG<N)>/(AM2-AM1> 

11 UH<l)=US<N)*Si8(l)*(lUG(NN-l)-UG(N))/DT 
AM1=(UU3(NN)+UB(1) )/2, 

SIB(2)=<1.-AMA*DT-XG<N)>/<AM2-AM1) 

. IF(A8S(SIB(2)-SIB(l))=TOL)30f30,12 
IP SI«( 1 )sS1B(2 ) 

GO TO 11 

3 n SSG ( NN ) = $G(N ) +S IB(2 ) * ( SSG < NN-1 > -SG ( N ) ) /DT 
IF < JCH-l)32,3i,3l 
31 A AG ( NN ) =UUG ( NN } 

O O TO 3 4 

3 P AAG(MN1=FXP< >*SS(i<NN)/2. ) 

34 AM1 = iJUG(NN)^AMG(NN) 

SIH(1)=(1..AM1*UT-XG(N)>/(AM2-AM1) 

13 U8(l)=UG(N)+SlB(l)*(UlG(NN-l)-UG(N.) ) / D T 
AB(1)=AG<N) + S*B(1 )*(AmG(NN- 1’)-AG<N))/DT 
AM1= C UllG(NN) + AAG( NN ) *UB( 1 ) +AB ( 1 ) ) / 2 . 
SIH(2)=(l,.AMi*DT-XG(N))/(AM2-AMl) 
lF(AHS(Slb(2)-SIBfl) )cT0L)15,l5.14 

14 SI«( 3 )rSIB(2 ) 

GO TO 13 

In SB(1)=SG(N)+SAB(2)* (SSG(NN-1 > -SG(N) )/DT 
XB<1>=XG(N)+A*2*SI3<2) 

DTP=DT-SI6(1 ) 

PB = 2.*A3(l)/(«-l. ) +UE ( 1 ) 

CALL AREA(XSd) .DAB) 

DDF = -UI)G(NN)*AaG( NN ) * D AE 
DDB=-UR(l)*ABll)*DAB 

PPG(NN)=PB^(DUE*DDB)*DTP/2.+( AAG(NN)+AB< 1) )*(SSG(NN)-SB(1) 5/2. 
IF! JCH-l)35,3o, 36 

3=5 IMJG(NN)=PPG(N'N)-2.*AAG(NN)/(G-1. ) 

GO TO 37 

3 6 t? UR ( NN )=(G-1.>*(PPG(NN>/(G+1.)) 

37 CONTINUE 

T F < ABS(UUG<NN)-UGUESS >-T Ol >15*16, 10 

16 IF (JCH-1 ) 2 3 # 1 * * 19 
27 XGR< NNS=XB(1 ) 

XGG C NN ) =G . 

IF ( UUG t NN ) .AAU ( NN ) )24 , 17 , 17 

17 JC^=1 

UUG ( NISI ) =UUG ( N'V-1 ) 

UGUESSr-JUG ( NN ) 

GO TO 2 n 

19 XGR(NN) = Xb(l ) 

X G 0 ( N N ) = 1 . 


24 OOG(NN)=2.*AA^(NN)/(Gcl, ) & U U G < N N ) 

RETURN 

END 
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3200 FORTRAN ( 2 . 1 . 0 ) / ( RTS ) 


SUBROUTINE: GEUM(ME>XE,DE»EE»C) 

DIMENSION xEf^b) # DE ( 25 ) , EE ( 25 ) * C ( 25 , 3 ) 

M 1 = M E - j 
DO 31 1=2, Ml 
J= I -1 

C(I,2)=(DE(I)-DE(J))/(XE(I)-XE(J>) 

C( 1 , 1 ) = UE ( J )-CE( I )*C( 1,2) 

3i C< ! ,3) = ( (DE( I*1 >-DE( I ) )/(XE( 1*1 ).X£( I ) )-C( I .2) >/<4.*EE( I > > 
RETURN 
FND 


SUBROUTINE ARgA ( X , DA ) 

CQMMlIN IL..IU > TqL .-I-OL-Q 

COMMON XE<25;,DEC2'?),ePS(25),YC(25*3>,M,XT ,XTR,XD 
IXLX-XTL ) 26. ^6 > 23 . 

26 rA=0, 

-C.=!1EX14 

r d = o . 

OH— m ±2 ; 

27 CONTINUE 

__ IFO.-1. )29.28.28 _ 

28 C A=0 . 

CjUIEIHI 

r d= o . 

GO TQ . 4 .2 

29 M1=M+1 

CQ. .30 1 = 2_,J11 

IF < X-XF < 1 ) >33,33,30 

30.. CONT INUE 

32 nn=(DE( I-1)-Dg( I ) )/(XE( I -1 ) -XE ( I ) ) 

H=JJE ( I -1 ) + DD» ( XrJl£ll-.llT. 

CO TO 46 

vB = XE( I >-EPSU) 

IfLl XA-.X1.4J. , .3 2^42 

41 TF(XB-X>44,3£,32 

—42. K = 1-1 

CO TO 45 

— 4.4- (C=T 

45 S I =X- ( XE(K ) -PpS( K ) ) 

fL=.Y Cl K.^LL*XCiK- , 2 ) *.SJ + Y.C.UU 3) pS-1*.*2. 

ri) = YC<K,2)*2**YCfK.3)*SI 

-.46 rA = 2..*flilZfl 

47 CONTINUE 

ELEJ.UKN. 

fnd 
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3200 FORTRAN < 2 . 1 . 0 > / < R TS > 

SUBROUTINE CRPSS( IQ.KO, JO, JJQ,NEWQ) 

COMMON G,DT , TOL . TOLO 

COMMON Xt(2“>)»DE(25).EPS(25),YC(25.3),H,XTL,XTR»XD 

COMMON XXG< 350 ) , PPG (350 ) » QQG( 350 ) . AAG( 350 ) ,UUG( 350) , SSG( 350 ) 

COMMON XG< 35 0> , PG( 350 ) , QG(350 ) , AG (350 ) , IJG< 350 ) , SG( 35 0 ) 

COMMON Z3 ( 350 > , ?2(35C ) , XGB( 35 0 ) , XGQ(350 ) 

CO M MON ?Ml,PWi,PPMl*FPWl 

COMMON PkFt , A1 ,U1 , PI, 01 , Si, A2, U2,P2,02,S?, A3,U3,P3,Q3. S3 
COMMON NNQ< 10 > , QQW< 10 ) , QGM ( 10 > » NO ( 10 > . QW <10 > , QM<10 > . K-NQT, NOT 
COMMON MmCS ( 1 o ) , MMCST , MCS ( 3. 0 ) , MCST 
COMMON XFHOCKilQ) 

N E W 0 = 0 

IF( JJQ-2)12.1«.10 
in K1 =NNU ( JJQ-1 ) 

IF(AriS(XXG(Kl)-XXG(IC>)-. 004)11. 11. 12 
11 10=10-1 
GO TO A 
1? continue 

no 1 1=1, NQ1 

K? = N()( 1 ) 

IF(A8S(XG(KC)-XG(K2))fc. 0 04)2, 2,1 
1 CONTINUE 
GO lu 3 
2 1.0=10-1 

x y t-i ( i . Q ) = y x G ( IN) 

PPG ( LQ ) =PPG ( TP ) 

OQi, ( I.Q ) =OOG (IP) 

A A G ( L Q ) = A A G ( I U ) 

UljG ( LQ ) =UUG ( I P ! 

SSG ( till ) =SSG ( ! p ) 

XGM(lG) = XGB( I p > 

X ! j 0 ( L Q ) = X G 0 ( IP) 
t Q = lO 
GO TO 4 

3 DQ= ( QQG ( 1 Q ) - Q P G ( IQ-1) ) / A A G ( IQ-1) 

0 0 M < J J Q ) = 1 , 1 

CALL SHOCK (CO* 03M( JJG ) , UR AT . AR A T, HR A.T , DS ) 

CQ'.- ( J JO ) = U U G ( l O-i ) - A AG ( IG»1 ) *QOM ( JJQ > 

X X' G < I Q > = X X G ( I P - 1 ) 

NNO < JJO ) = I 0-1 
M p *v Q = 1 

V SHOCK ( JJQ ) = X*G (IQ) 

JJO= JJO+1 
4 CO Ml INI IE 

PET URN 
FNO 
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3200 FORTRAN (2.1.0)/ ( R r S ) 

S'JKROUTlNfc G CROSSCUT. I 0 » KQ , JO , J JO . JCS , J JCS , I1Q.KKQ.N) 

COMMON S.OT.TOL.TOLQ 

COMMON XM 25 > » DE(25) , EPS ( 25) , YC ( 25 . 3 ) , M , XTL , X TR , X D 

COMMON XXM35U) , PPG (350 ) ,QOG( 350 > ,AAGC 35(1 ) ,UUG( 350 ) ,SSG< 350 ) 

COMMON x G ( 3 5 0 > , P G ( 3 5 0 ) >QG(350 ) » A R ( 3 5 0 ) » U G ( 3 5 0 ) .SG (350 ) 

COMMON 21 (350) , Z2C350 ) . XGb( 350 ) , XGQ(350 ) 

COMMON “M1.PWA.PPM1.FPW1 

COMMON P2P1 .Al.Ul.Pl.Ol.Sl, A2.U2.P?,02.S2.A3.U3,P3.03,S3 
COMMON NNCldO ( . QQW(10 ) . QGMdO ) .NGK10 ) .GW(10 ) , QM(10 ) , NNGT. NOT 
COMMON MMCS(1U).MMCST,MCS(10),MCST 
COMMON X SHOCK dO ) 

CO M MON XCONT A“ T ( 10 > ,L CONTACT (10 > 

K 1 = JiJ -V 
K2= JO-1 
K 3 = N Q ( K 1 ) 

K4=K3+1 

K5=NQ(K2) 

K6=K5+1 

AMI = r 00 w( JJC-1 )+OW( JGel ) )/2 . 

AM?=(OOW( JJC-O+0W< JGe2) )/2. 

TIN=( X SHOCK ( jwo-l )-XSH0CK(JJ0-2))/(AM2-AMl)+T 
XlN=XSHDCK(oJW-l ) + A M 1 +■ ( T I N - T ) 

SS4=SG(K6) 

PQ/i = QG(<6) 
lllj5 = UG ( KA ) 

A MG!l£SS = 005 

00= ( UU5-UG ( K3 ) ) / AG ( K3 ) 

0M5=(-0J*(G + 1 • )/2.+SCRT< UJUMG+l. )/2. )**2 + 4. ) )/2. 

C A I L GhOCKO ( O 05 , DU5 , AR AT . PR AT , DS . DO ) 

AA9=ARAT* AG( K«J ) 

S S 5 = S G ( X 3 ) + C S 

AA4 = AA5*FXP< (o-l. )*(SM-SS5)/2. ) 
ilil4 = ?.*AA4/(6"l. ) - QQ4 
1 1 IJ 5 = 1 1 U 4 

IF ( A8S< 005. L GO ESS) -TCL ) 7, 7, 6 
7 nw5=OG(X3)-AG(K3)*0M5 
XS = X f N+0W5+ ( T-T I N ) 

VCS=X 1 N+UU5* ( I -TIN) 

DO 8 1=1,10 

I F ( X X G ( i )-XSHUCK( JJO-1) ) 8 , 9 , 9 
A CONTI Nile 
1 0= 10+1 
GO rn j u 
9 10=1 

in V X t. ( l Q ) = X S 

PPG ( I0)=P6<«3> 

P0G< I 0 ) = 0 G ( K 3 ) 

AAG( I Q ) = AG( K3 > 

HUG (IQ) = U G ( X 3 ) 

SSU< IG)=SG(<3> 

XG«( I Q ) =0 . 

X60 < I Q ) = 0 . 

JJO= JJO-2 
X SHOCK ( oJ0) = X“ 

NNO( JJO)=IO 
00 w ( J J 0 ) = 0 W 5 
OQM ( JJO ) =QM5 
JJO= JJ0+1 
10= I o+l 
XXG( I Q ) =XS 
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A A G ( ! (J ) = A A 5 


uunt IQ i=HU5 

P P G ( I U ) = 2 . ★ 4 A o ( j Q ) / ( 0 k 1 . ) + U 0 G ( IQ) 
OQi. ( I Q >=',».★ A A* ( IO)/(Gr.l. )=UJG( IQ) 
^ S L- < IQ ) =SS5 
V(JH( I Q ) = 0 . 

V G 0 ( I Q ) = f) . 

10= I 0+1 

MHOS ' JJCS) = f 0 
yy(- ( I 0 ) =vGS 
ppr, ( IQ) = ppg< iw-i) 
oqo( i q ) = r: o fi ( iu-D 
A AP( T 0 > - A A G ( IvJ-1) 

DUG ( I Q ) = I ' U G ( I J - 1 ) 

S5«( 1 Q ) =SSG ( W-l) 

VG M ( I Q ) =n . 
y G o t I a ) = u . 
yCONTACI ( JJCSJ=XCS 
l.lCONT APT ( J JCS > =UU4 
JJCS=JJGS+1 
J U= 1 Q + 1 

y x g ( IQ) = ycs 

A A G ( I Q ) = A A 4 
DOG ( ! Q ) = D U 4 

PPM IQ) =^.*4 AM IQ)/(G?1. )+UiJG( IQ) 
OQ(i(IQ)=2.*AAu(IQ)/(Gsl.)sUiJG(IQ) 
SSG( IQ ) =SS4 
XGH f IQ ) =(i . 

X G 0 ( IQ) = n. 

T I o = I Q 
K K 0 = * Q 
P E T ij P |\j 
PHD 
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3200 FORTRAN ( 2 . 1 . 0 > / ( R TS > 


SUHRQUTIN& W 0 V E ( N N » N * 1 PGO 3 

common h.dt.tul. tolq 

COMMON X f ( 2 9 ) « 0 1 ( 2 5 ) . E:PS(25),YC(?5*3) »M,XTL,XTR«XD 

COMMON XXti(i5U > f PR6(J£iO ) . OOQ(35l)> t AAG(3S»0 > .UUG( 35U ) ,SSG( 350 > 

COMMON XG<350 > . PGC350 > ,QG< 350 > » AG<350 > »UG<350 > » SG< 350 > 

COMMON 1 1 ( 350 > , '/.?( 350 ) , XGB (350 ) , XG'K 350 > 

COMMON 3 Ml , PWA , PPM1 , FPW1 

COMMON p?P 1, Ai.ijj ,Pl.Gl,Sl,A?»U2,P?,02»S2»A3.U3,P3»Q3iS3 
COMMON NNO< 10 t , QQW( 10 5 , ROM (10 ) * NO U 0 ) , GW (10 ) , QM(10 ) , NNQT , NOT 
COMMON MMCSdU ) , MMCST , MCS( 10 ) , MGS I' 

no 14 I = 1 » NN 

X G < I ) = X X G ( I ) 

■ PG ( I ) = P P G ( I ) 

QG ( I ) = OQG ( I ) 

A G ( I ) = A A G ( I ) 

! ! G ( I ) =IHJG( I ) 

14 R G ( l ) =SSG( I ) 

N = N N 

NQ >' = N N C T 
no 1.5 1 = 1 # N G T 
NO( ( )=W.MO( I ) 

r>w 1 1 ) = paw ( I ) 

15 PM ( I ) =QUM ( I ) 

NQT1=N0T+1 

!'0 ih I=NUTl,iO 

1 M O ( [) = •'! 

M C S I = H M G S T 
HO 1 / I = 1 , M C S 1 
17 MCS ( I) =MMCS ( I > 

M 1 = M C S T + 1 

no 2 0 i = f ; 1 1 1 o 

2 n M C S ( I ) r ii 

I F < JOG 0)19, 18*19 
1 A PMlsPHMl 
PW1 . sPPWl 
19 CONTINljt 
RETURN 
END 
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3200 FORTRAN ( 2 . J . 0 ) / ( RTS ) 

SUBROUTINE CS^POSSQ( T, IQ.KO. JO, JJQ. JCS. JJCS . I IO.KKQ.N ) 

COMMON i,DT,TRL, TOLQ 

COMMON XE(25)'DE<25),EPS(?5),YC(?5>3),M,XTL,XTR,XD 

COMMON XXG(35U ), PPG (3 5.0 > ,QQG< 350) »AAG <350 > »UUG(35U > ,SSG( 350 > 

COMMON XG ( 35 0 > , p G(350 ) , QG< 350 ) , AG (350 ) . UG(350 ) > SG( 350 ) 

COMMON Z1 (350 > , Z2(350 ) , XGB<350 ) » XGO(350 ) 

COMMON Phi ,PWi , PPM1 ,FRW1 

COMMON P?P1,AX,U1,P1,01,S1,A?,U2,P2,02.S2,A3,U3,P3,Q3,S3 
COMMON NNU< 10 > , QQW(1C ) , QOM( 10 ) » NQ( 10 ) , QW<10 ) , QM<10 ) , KNOT, NOT 
COMMON mmCS( 1U),MMCST,MCS(10),MCST 
COMMON XSHOCKtlO) 

COMMON XCONTAPT (10 ' .LCONTACT (10 ) 

K1=NG(J0-1) 

K 2 = M C S ( JCS-1 ) *1 
K3=MCS< JCS-1 ) 

K 4 = * 3 + 1 
K5 = NO< JO-1 ) 

K 6~K5+i 

OX H T CS = ( U C 0 NT A C T ( JJCShl ) + U G ( K 2 ) )/2. 

PXOTS = Ori( JO-l > 

TIN=T-DT+(XG( *1 ) -XG(K2) ) / ( DXDTCS-DXDTS) 

X I N = XG ( M )*CXMTS*( T INsT + DT) 

IF(TIN-T)V,fi,0 
A T I N = T 
Q SS4=SG(*6) 

OQ4=OG<*6) 

HU5=UG( <6 ) 

If) MGU£SS=UU5 

DU= ( UU5-UG ( K 3 > ) / A G ( K 3 ) 

QM5 = ( - O J * (G + l • ) /2 . +SCRT ( ( DU* ( G + l . )/2. )**2 + 4. ) )/2. 

CAl L SH0CK0(0rl5,DU5. AR AT , PR AT , D.S . DQ ) 

AaS=ARAT*AG(KP) 

SS5=SG(*3)+CS 

AA« = AA5*FXP( (»»-l. )*(SS4-SS5)/2. ) 

UU4 = 2 . * A A 4 / ( G"1 . ) -QQ4 
UU6=UU 4 

I F ( ARS( UU5-U GUESS ) -TCL ) 11, 11, 10 
11 OW5=UG(<3).AG(K3)+QM5 

PRPLOLO=?.*G*«M( JQ-1)**2/(G*1. )-(G-l.)/(G+l.) 

IF(PHAT-PRP|_0l-n)2,2»l 
1 WR I Th ( 61 , 1 0 0 ) 

ion FORmaT(3X,33HpONTaCT SURFACE CROSSES A 0 SHOCK/) 

0MN6W=$QRT ( (G + l. )*( (AA4/AG<K6) )**(2.*G/(6”1. ) ) + ( G - 1 . ) / ( G + l . ) ) / ( 2 . * 
1G) ) 

WR I TE ( 61 , 101 )UMNEW 

10i FORMAT(3X.3lHlHE REFLECTED RAVE IS A P SHOCK, 2X»5HMACH=, El5.fi ) 

WR I TE ( 61 , 102 >MRAT 

10? FORMAT (3X.-33HHR6SS RATIO OF TRANSMITTED SHOCK = , El? . 8 ) 

W R i T t ( 61 ,103 )RRPLOLD 

103 FORMAT ( 3X,30HHRfcSS RATIO OF INCIDENT SHOCK =, E 15 , 8 ) 

WRI TE ( 61 , 104 ) 

104 FORMAT* 3X.56H I HE PROGRAM TREATS THE P SHOCK AS VERY WEAK AND DROPS 
1 IT) 

WR I TE ( 61 . 105 ) I IN 

105 FORMAT ( 3 X » 6 9 H i F THE F SHOCK IS NOT WEAK, ALL RESULTS ARE NOT CORRE 
1CT LATFR THAN T = ,.E15-8) 

WR I TE ( 61 , 11 0 ) 

110 FORMAT ( 1H1 ) 

? CONTINUE 

XS=X IN + Ok'5* ( T-T IN) 


XCS = X i N+UU5* ( I -T IN ) 

DO 12 1=1 , IG 

I F ( X X G ( I )-XSHUCK( JJQ-1) ) 12. 13, 13 
1? CONTINUE 
I Q = I Q + 1 
GO TO 14 

13 IO=I 

14 XXG ( I Q ) =XS 
PPG ( 1 Q ) =PG ( K3 ' 

OQG ( IQ>=GG(K3> 

AAG(IQ)=AG(K3) 

(JUG ( IQ)=UG(K3> 

SSG ( JQ)=SG(K3> 

X GN ( IQ’)=0. 

XGOf IQ)=0 . 

JJ0= JJO-l 
XSHOCK ( JJQ ) =X3 
MNO( JJ0)=IQ 
OOw ( JJO ) = QW5 
PQM( JJO)=QM5 
JJ0 = jjo + i 

I Q = 1 Q + 1 
X X G ( I Q ) = X S 
A AG ( I Q ) = A A ? 

IMJG ( I 0 ) =UU5 

PPG( I Q ) = 2 . * A A G ( JQ)/(G»1. ) + U U G ( IQ) 
OQG( IO)=?,*AA*J( IQ) / (Gal. ) bUUG (IQ) 
GSG( I Q ) = SS5 
XG8 ( IQ)=0. 

XGQ ( I Q ) = 0 , 

I Q= I 0 + 1 

JJCS = J.JCS-1 

MMCS ( J JCS ) = I 0 

X x G ( I Q ) = X C S 

PPG ( I Q ) = P P G ( IQ-1) 

OQG ( I Q ) =QUG ( IQ-1) 

A A G ( I Q ) = A AG ( IQ-1) 

UUG ( I Q ) =UUG (IQ-1) 

SSG ( I Q 1 = SSG ( IQ-1) 

XGH ( IQ ) =0 . 

X G 0 ( I U ) = 0 . 

XCONTACT (JJCS^=XCS 
UCONTACT ( JJCS > =UU4 
JJCS=JJCS+1 
10=10+1 
XXG ( I U ) =xcs 

A A G { I Q ) = A A 4 
IMJG C I Q 1 =UU4 

PPG( IQ)=?.*AAQ( lO)/(Gs.l . )+UUG( IQ) 
Q 0 G ( I Q ) = 2 . * A A u ( I 0 ) / ( G s 1 . ) a U iJG ( IQ) 
SSO ( I Q ) =SS4 
XGP( IQ ) =0 . 

XGO( I Q ) = (1 . 

I 10= IQ 

kko=kq 

RETURN 

END 
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3200 FORTRAN ( 2 . 1 , 0 ) / ( RTS ) 


SUBROUTINE SHUCKQtOM.UU, APA T , PRAT , US . DQ ) 

COMMON G 

nu=2.*(l. -QM**2) /( (G+l. )*QM) 

ARATsSORT((2. + (G-1.)*0M**2)*(2.*R* ( 3M**2-(G-1.)))/((G + 1.)*QM) 
PRAT=2.*G*OM* w 2/(G+l. ) - ( Gal . ) / ( G + l . ) 

AG1 =2 . *5*(JM**2/ (6+1 . (G-l . ) / ( G + l . ) 

A G 2 = ( (1 . + ( G - 1 • )*OM**2/2. )/( < G + l . > *GM**2/2 . ) )**G 
DS=AL0G(AG1+Au2)/(G*(C-1.)> 

DQ=2.*<ARAT-1. )/(G-l. )+DU 

RETURN 

END 


320 0 FORTRAN ( 2 . 1 . 0 ) / ( R TS ) 

SUBROUTINE GDiS( JO, XC 1S|_. XD I SR > 

COMMON G,DT.TUL.TOLQ 

COMMON XE(?5)»DE(25),EPS(25),YC(25.3>,M,XTL.XTR,XD 

COMMON XXG(35U ) .PPG(3J)0 ) ,QQG(350) . AAG(35H ) .UUG(350 ) *SSG(350 ) 

COMMON XG< 350 ) , PG( 350 ) ,QG< 350 > . AG< 350 > »UG< 350 ) »SG t 350 ) 

COMMON ZK350 J , Z? ( 350 ) , XGE < 35 0 ) .XGQI350) 

COMMON ?Ml,PWi,PPMl,FPWl 

COMMON P2P1 . AX, ui, PI. 01, SI, A2 . U2,P2,Q2«S2» A3,U3,P3, Q3» S3 
COMMON N N Q ( 1 0 > , QOW < 10 5 . QQM ( 1 0 ) , NO ( 1 0 ) , QW ( 1 0 ) , QM { 1 0 ) , NNQT , NOT 
IF ( NOT ) 1 , 1 , 2 

1 XDISL=0. 

XOISR=l.l 
GO TO 9 

2 IF(NGT-JC!)5,3»6 

3 I F ( JO-1 ) 4 , 4 . 8 

4 K1=NQ( JO) 

xnisi = n . 

XOISR=XG(Kl ) 

GO TO 9 

5 K2=NQ( JO-1 ) 

' XDISL=XG(K2) 

XD I SR=1 . 1 

GO TO 9 

IF ( JO- 1 ) 7 » 7 * 8 

7 Kl=NQ(ja) 

XD I SL = 0 . 

XDISR=XG(K1 ) 

GO TO 9 

8 Kl=NO ( JQ ) 

K2=NQ( JQ-1) 

XD1SL=XG(K?) 

XDISR=XG(Kl) 

9 CONTINUE 
RETURN 
END 
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3200 FORTRAN ( 2 . 1 . 0 ) / ( RTS ) 


SUBROUTINE X D l SCON ( JG i JCS , X D I SR , NQOCS ) 

COMMON G.DT.TUL.TOLO 

COMMON XE(25 ) . DE(25) . EPS ( 25 ) , YC ( 25 . 3 ) ,M,XTl,XTR,XD 

COMMON XXGt 35 U ) , PPG ( 350 ) , OOG ( 35 0 ) , A AG ( 350 ) ,UUG (35 0 ) , SSG ( 350 ) 

COMMON XG(35 0 > .PGC350 > , QG(350 ) , AG (350 ) .UG(350 ) , SGi 350 ) 

COMMON Z1C350>, 22(350), XGB ( 350 ) . XGQ ( 350 ) 

COMMON °Ml,PWl.PPMl,PPWl 

COMMON »?P1.A1,U1,P1.Q1.S1,A2,U2,P2,Q2.S2.A3,U3,P3,Q3,S3 
COMMON NNQ ( 1 0 ) , QOW(10 ) ,OOM( 10 ) ,N(3(10)»QW(10),QM(10)»NNQT,NQT 
COMMON MMCS(IU) .MMCST ,MCS(10) ,MCST 
IF(MCST)3,3,1 
3. IF (MCST-JCS)3>2#2 

2 KlsMCS ( JCS > 

XRCS=XG(K1> 

GO TO 4 

3 XRCS=l.i 

4 IF ( NOT - JQ ) 6 > 6 • 5 

5 K1 = NQ ( JQ + 1 ) 

XRQS=XG(K1> 

GO TO 7 

6 XRQS = 1 . 1 

7 IF(XRCS-XROS)10.9,8 

8 XDISR=XROS 
NQOCS=l 

GO TO 11 

9 XDISR=XRCS 
NQOCS=0 

GO TO 11 

10 XDISR=XRCS 
NOOCS=2 

11 CONTINUE 
RETURN 
END . 


3200 FORTRAN < 2 . 1 . 0 ) / ( RTS ) 

SUBROUTINE SHOCK ( DEI , SM , UR , AR , PR, DS ) 

COMMON G.DT.TOL 
Sl=-2. 

52 = ( G*1 . ) * ( DEU+2 . / ( G- 1 . ) ) 

53 = 2 . 

R4=8.*G/(G.l . ) 

S5=16 , *G/ ( G-l • )**2-4. 

S6 = -8./(G-l . ) 

T1=S1**2-S4 

T2=2.*S1*S2 

T3=2.*S1*S3+SK**2-S5 

T4=2.*S2*S3 

T5=S3**2-S6 

1 SMG=SM 

F=T1*SM»*4*T2*SM**3+T3*SM**2+T4*SM+T5 

FPr4.*Tl*SM**P+3.*T2*SM«*2+2.*T3*SM+T4 

RM=SMG-F/FP 

IF(ABS(SM-SMGI-TOL)2.2»l 

2 CONTINUE 

UR=ABS(2.»(1.-SM**2)/((G*1,)*SMM 

AR = SQRT( (2.* (0-1. )*SM**2)*<2. *G*SM**2- (G-l. ) ) )/< <6+1. )*SM) 
rR = 2.*G*SM**2/(G + l . )- ( G-l , )/ ( G + l . > 

AG1 = 2.*G*SM**4/(G + 1. ) » ( G ■ 1 . ) / ( G + l . ) 

AG2=( ( 1 . ♦ ( G- 1 • )*SM**2/2. )/( (G+l. )+SM**2/2. ) )**G 
DSsAL0G(AGl*AO2)/(G*((i-l. )) 

RETURN 

END 


3200 FORTRAN ( 2 . 1 . 0 ) ✓ ( RTS > 


SUBROUTINE PHUCK< IQ,N,NN> 

COMMON G,DT.TUL,TOLQ 

COMMON XE(?5 ) / DEC 25) , EPS (25) , YC ( 25 , 3 ) , H , XTL . X TR , XD 

COMMON XXG(35U ) .PP6( 3 50 ) , COG (350 ) » AAR (350 > .UUG< 33d ) , SSG( 350 ) 

COMMON XG( 350 * .PG(350 ) , QG( 350 ) » AG ( 350 > > UG ( 35 0 ) » SG ( 350 ) 

COMMON ZK350 > . Z2< 350 > , XGB< 350 ) . XGQC 35 0 > 

COMMON »Mi,pwl # PPMl f FPWl 

COMMON P2P1, Al.Ul, PI, 01. SI. A?,U2.P?,Q2>S2, A3.U3,P3,Q3,S3 
DIMENSION ClU).XB(l).Ufi(l),AB(l).SB(l) 

PWl=PMl 
SP=PM1 
in SGG=SP 

UAsUl+?.*Al*(PP**2-l.)/(<G+l.)*SP) 

AA=SQRT{(2.+<^-l.)*SF**2)*(2.*G*SP**2-(G«l.)))/((G*l.)*SP) 

PPW1=U1-A1*SP 

AG1 = 2.*G*SP***/(G«-1. >s<G-l. >/<G*l. > 

A G2= ( ( 1 . ♦ ( G- 1 • )*SP**2/2. ) / ( ( G + l . )*SP**2/2» > >**G 
SAsSl^ALnG(AGX*AG2)/(G*(G«l.)) 

XA = XG(N)«-DT* (PWl + PPWl ) /2 . 

CALL AREA(XA.UAA) 

Cl(l )=XA-DT*(^A + AA)/2 , 

CALL IWTFRf2,!M,i,zl#XG*Cl*XR) 

CALL ARE A ( XB ( 1 ) , D AB ) 

CALL lNTEH(2#N,i,xG#Lli»XB 4 UB) 

CALL lNTfcR(2,i'i,l,XG#AG,XB J A9> 

CALL INTER (2 , * , 1 , XG, SCI, XB, SB) 

PBe2.*A9(l)/(V»-l. ) + UE ( 1 ) 

PAsPR-M - UA*AA*DAA-UB( 1 > *AB< 1 >*DAB >*DT/2 . *< AA + AB( 1 ) )*< SA-SB ( 1) > /2 . 
DEL = < P A - PI ) / A 1 

CALL SHOCK ( CEL, SP , URA T , A P A T , PR, DS) 

IF( ABS(SP-SGG>-TOL)ll j 11* 10 
11 CONTINUE 
NN= I Q + l 
XXG(NN)=XA 
PPG ( NN ) =P A 
DUG ( NN ) =11 A 
A AG(NN ) rAA 

QQG(NN)=?.*AA/(G-1. ) - U A 
SSG(NN)=SA 
X GB ( NN ) = X B ( 1 ) 

V G 0 ( N N ) = 0 . 

PPMi=SP 

RETURN 

FND 


3200 FORTRAN < 2 . 1 . 0 > / < R TS > 

SUBROUTINE 1NIER<N2,N # N1.X»Y,XS,YS> 

DIMENSION X(3u0),v(300>.XS(l)tYS(l) 

IF(XS(1)-X(1> >1.1,2 
i J = 2 
GO TO 7 
P DO 3 1=1, N 

IF { X S ( 1 )-X(I ) >4,4,3 
7 CONTINUE 
J = N 

GO TO 7 
4 J=I 

IF(ABS(X(J)-X*J-1>>-. 000001)6, 6. 7 
ft YS ( 1 ) = ( YS ( J ) + Y S ( J-l ) ) /2 . 

GO TO 8 

7 YS(i) = Y(j-n^vs(i)-xtj-i))*(Y(j)-Y(j.i))/(x(j).xu-in 
a CONTINUE 
RETURN 
FND 

3200 FORTRAN DIAGNOSTIC RESULTS - FOR INTER 
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3200 FORTRAN ( 2 . 1 , 0 ) / ( RTS > 


SURROUT 1 NE PRINTOUT ( T , NN. IPGO.XTEST) 

P I MENS I ON XT(1).UT(1),AT<1).ST(1> 

COMMON u,L)T,TUL,TOLQ 

COMMON X6(?5)»DE(25) ,fcPS<25), YC(25,3),M,XTL.XTR,XD 

COMMON XXG ( 350 ) , ppG ( 350 ) , QQG ( 350 > , A AG ( 350 > . UUG < 350 ) , SSG ( 350 > 

COMMON XG(350>,PG(350).GG<350 ) , AG ( 35 0 ) . UG ( 350 ) »SG(350> 

COMMON Zl(350>, 22(350), XGb(350),XGQ (350) 

COMMON PMl,FWl,PPMl»FRWl 

COMMON P?P1,A1,U1,P1,01,S1, A2,U2,P2,Q2,S2, A3.U3,P3,Q3,S3 

common nngkio < ,oow(io ) , qqm(io ) ,no do ) , qw<io > ,qm(io ) ,nnqt, not 

COMMON MMCSdU ) ,MMCST.MCS(10 ) , MCST 

common xshock ' It) ) 

COMMON XCONTAUTdO) .ICON TACT (10 ) 

OP = G/(G-l . ) 

XT d ) =XTEST 
WR1TE(6 1,135)' 

115 FORMAT ( //3x»5«TlMe= , E15. 8, IX . 3H( . ) // ) 

WR! TE< 61,117 ) 

117 FnRMAT(5X,lPI,7X,3HX,14X,lHP,14X,lrtQ,l4X,lHA,14X,lHU,14X', 1 HS, 1 UX ,8 
IMP ORIGIN, 8X, OHO ORIGIN/) 

JJO = 3 

JJCS=1 

DO 8 1=1, NN 

WRITE (6 1.1 (13 H,XXG< I ),PPG( I ) ,G)QG( 1 ) , AAG( I ) ,UUG< I ) ,SSG ( I ) ,XGB( I ) , 

1 X G 0 ( 1 ) 

103 F ORM AT ( 3x , I3,t»El5.8) 

IF(NNQT)5,5,2 
2 I F ( J JO- NMQT ) 3 » 3 , 5 
7 IF( l-NMQ(JJC) >5.4,5 

4 WRIT 6(61. 102 )<J()M( JJO),OQW( JJQ) 

102 FORMAT ( 3X , 3 7HW SHOCK MACH N0. = ,E15.8,4X,17HQ SHOCK VELOC I TVs , E15 , 8 
1 ) 

JJU=JJO+l 

5 IF(MMCST)6,e,o 

A IF ( JJCS-MMCST ) 18, 15, e 
15 IF( I -MMCS ( J„Ci> ) )8,7,e 
7 WR]TE(61,B7) 

87 FORMAT f 3X . 15HP0NT ACT SURFACE) 

JJCS=JJCA*j 
A CONTINUE 

I F ( IPG 0)10, 9, AO 
9 WR I'l'R ( 61 , 135 ) MPM1 , ppw 1 

135 Fow MAT ( / 3 X , 2 0 1-1 d SHOCK MaCH NUMBERS , El5 . 8, /3X , 17HP SHCCK VELOCITYs, 
1615.8) 

10 CONTINUE 

IF ( XXGd ) -XT (X) ) 11 , 12 , 12 
11 CONTINUE 

CALL I N T bR ( 2 . r«N , 1 , X XG , UUG , X T , UT ) 

C At L INTER(2,'VN,1,XXG, AAG,XT,AT) 

C A I L lNl'ER(2,i'tN,l,XXG,SSG,XT,ST) 

AMACH = l!T ( 1 ) / A i (1 ) 

PR6S = ATd)**(C.*G/(G-l.))*EXP<-G*ST{l>) 

WRITE (61. 100 )XT( 1) , APACH.PRES 

10n F0PMAT(//3x,6MxTEST=,E15.S,2X,5HMACH=.El5.8,2X,6HRRESS=,El5.e) 

17 WRITE(ftl,110) 

11 n FORMAT ( lhl ) 

PETURN 

END 
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